J.

HES;

 

IPLHLI .

..mniu..flz...ui.

. 2 2.!

53?. y 4314.”.

. . . .1 $58.53,:
_ is

5:42...-‘1
I... II .4
$23!...
6.1». :1:
..

. I}: J .v. I
£40440}. molt.)
.2») .1 ..YE.

‘5...)91...‘ $5.):
365.155....

l‘alx‘:&
1.‘

t3. ..
..l . ‘sifi?
pita—‘15:.

fivdhffifi.

.75 1 3311.411}! EC:- ‘5
$35.35.!IP5EZ‘. .t-IE; .
. .Esoxi..¢u.vmm.”“i fit-.5“
9|.
. .19.
[-51.0155

 

..::$1... .. . . . . . , . . . . . . 933... .. .......n....ué.u....:...1
.. . . . , H . . _ . . _. _. . . . ..wfigmafia.

 

 

 

LIBRARY
Michigan State
“ ‘ University

9?C_\A

s.__/'

 

 

 

This is to certify that the
dissertation entitled

MECHANICAL MODELING AND SIMULATION OF POROUS
POLYMER NETWORKS: LOAD INDUCED LOSS OF
SATURATION IN ISOTROPIC ELASTOMERS AND
PRESSURE DRIVEN SEEPAGE IN DIRECTIONALLY
REINFORCED ELASTOMERS

presented by

HUA DENG

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Mechanical Engineering

 

 

WW...

- Major rof firs Signature
:Dp adv/Aw / 6:. 200?

Date

MSU is an Affirmative Action/Equal Opportunity Employer

 

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES retum on or beforedate due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 KIProj/Acc&Pres/CIRCIDateDue.indd

 

MECHANICAL MODELING AND SIMULATION OF POROUS POLYMER
NETWORKS: LOAD INDUCED LOSS OF SATURATION IN ISOTROPIC
ELASTOMERS AND PRESSURE DRIVEN SEEPAGE IN DIRECTION ALLY
REINFORCED ELASTOMERS

By
HUA DENG

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Mechanical Engineering

2009

ABSTRACT

MECHANICAL MODELING AND SIMULATION OF POROUS POLYMER
NETWORKS: LOAD INDUCED LOSS OF SATURATION IN ISOTROPIC
ELASTOMERS AND PRESSURE DRIVEN SEEPAGE IN DIRECTIONALLY
REINFORCED ELASTOMERS

By
HUA DENG

Elastomeric gels are high molecular weight crosslinked polymer networks immersed
in a low molecular weight liquid medium. In the liquid environment, they could un-
dergo a large deformation associated with swelling or shrinking in response to environ-
mental stimuli, such as change in temperature, chemistry of the liquid bath, and light
exposure. This valuable property makes them useful in a wide range of applications in
drug delivery, surgical dressings, artificial tissue, and control material in engineering,
which motivate the desire to better understand their underlying mechanical response.

In gels the polymer and liquid components mix in definite proportions as de-
termined primarily by entropic and enthalpic effects. Mechanical loading can also
alter the mixture proportions by absorbing or driving out the liquid. Gel swelling in
the absence of mechanical loading is often described by a generalized Flory-Huggins
equation, which accounts for the effects between such a treatment and the broader
hyperelstic theory which accounts for the effect of mechanical loading. In this study
we consider loadings that can lead to both fluid gain (swelling increase) and fluid loss
(swelling reduction). For loadings that give fluid gain, we then consider a situation
in which the amount of available fluid is limited. In this case, increased loading may
reach a point at which no additional fluid is available for uptake into the gel system.
This results in a transition of the gel from a state of liquid saturation to a state
in which it is no longer saturated. This transition is first considered in the context

of homogeneous deformation where an appropriate hyperelastic analysis shows that

the transition from saturation to nonsaturation gives rise to an abrupt mechanical
stiffening. Then two kinds of inhomogeneous deformation problems are investigated,
including everting an axially loaded tube and twisting a hollow tube that originally
swells freely in the liquid bath. Various boundary displacements and traction condi-
tions are applied so as to study how these alter the original fluid distribution. It is
found that certain boundary conditions generate an overall volume increase after free
swelling, which results in a stiffer mechanical response after loss of saturation.
These static problems describe equilibrium situations in which both the fluid com-
ponent and the polymer matrix component of the system are at rest. However, a more
complicate phenomenon - which attracts abundant research interests - fluid diffusion
through polymer networks requires further study of the relative motion between the
fluid and the polymer network. A mixture theory is then invoked to specifically deal
with separate mechanical balance principles for each component. Pressure driven
fluid seepage problems for both isotropic and anisotropic (fiber reinforced) gels are

discussed based on this framework.

To my family

iv

ACKNOWLEDGMENTS

First I would like to take this chance to thank my major advisor, Dr. Thomas
Pence, for his guidance, instruction and support through the course of my doctoral
study and research. His expertise in solid mechanics and related areas makes me feel
easy to handle any difficulty during my 4 years Ph. D. career. He is always willing
to help me out as much as he could not only in my academic research but also in
my personal affairs. I would also like to thank my advisor Dr. Neil Wright for his
support, in particular valuable suggestions on the details of my dissertation format.

Second, I want to thank Professor Ronald Averill, Professor Xiaobo Tan, and
Professor Seungik Baek for serving on advisory committee. Their insightful comments
on my research are greatly appreciated. Professor Seungik Baek’s experience in the
similar research field has helped me a lot in completing the theoretical work of this
study.

Third, I want to thank my family for giving me the chance to experience this tran-
sient, but wonderful journey of life in this world, even though not much achievement
has been made yet.

Finally, I also would like to thank friends at Michigan State who have offered me
generous help: Dr. Li Sun, Guokai Zeng, Dr. Yang Fang, Fang Xie, Dr. Tengfei Luo
and many others. I want to express my special thank to Ms Manshan Li for her love

and support in my life.

Table of Contents

List of Figures ................................. viii

1 Introduction ............................... 1

1.1 Background .............................. 2

1.2 Objectives and Scope ......................... 11

2 Modeling of the Swelling Behavior of Elastomeric Gels ...... 15

2.1 The F lory-Huggins Equation for the Determination of Fme Swelling 15

2.2 Hyperelastic Constitutive Theory .................. 22

3 Homogeneous Deformation ...................... 29

3.1 Equibiaxial Stress ........................... 30

3.2 Uniaxial Stress ............................ 34

3.3 Experimental Validation ....................... 39

3.4 Equitriaxial Stress .......................... 41

3.5 Some Energy Considerations ..................... 47

3.6 Summary ............................... 52

4 Eversion of a Swollen Cylindrical Tube ............... 55
4.1 The Inhomogeneous Deformation of an Everted Tube under Axial

Load .................................. 56

4.2 Numerical Results for a Saturated System ............. 61

4.3 The Everted Cylinder that is Not Saturated ............ 68

4.4 Summary ............................... 79

5 Twisting of a Cylindrical Tube .................... 80

5.1 Twisting of a Saturated Tube .................... 81

5.2 Twist at Constant Overall Volume ................. 86

5.3 Twist Induced Fluid Desorption ................... 89

5.4 Twist Induced Fluid Absorption ................... 96

5.5 Fluid Uptake Leading to Loss of Saturation ............ 107

5.6 Summary ............................... 115

6 Seepage Through a Fiber Reinforced Hyperelastic Layer ..... 118

6.1 A Biphasic Mixture Theory ..................... 119

6.2 Free Swelling of a Fiber Reinforced Hyperelastic Layer ...... 126

6.3 Constitutive Theory ......................... 130

vi

6.4 Steady State Diffusion through a Fiber Reinforced Layer ..... 132

6.5 General Formulation for Non Steady State Diffusion ........ 150

6.6 Summary ............................... 172
7 Conclusion ................................ 174
APPENDIX .................................. 177
A Nondimensional Form of Governing Equations of Steady State

Diffusion Problem ............................ 178
B Functions F1(g’) and F2(g’) ....................... 179
C Function F3(f',f”,f’”,g’,g”,g’”) .................... 180
Bibliography .................................. 182

vii

2.1

2.2

2.3

3.1

3.2

3.3

3.4

3.5

3.6

List of Figures

The constitutive function h(J) in (2.8) for various X. Vaues of J giving
h(J) = 0 model gel dilation on the basis of mixing entropy and mixing

enthalpy, but do not account for the elastic effect of polymer crosslinking. 19

The constitutive expression h(J) + pJ-1/3 for various X when p =
0.01M. The new term pJ"1/3 describes the effect of crosslinks. Values
of J giving h(J) + uJT1/3 = 0 describe free swelling. .........

The free swelling stretch ratio C as given by (2.14) and (2.15) is a
function of M, X and E. The dependence of C on g is shown for various
M and various X .............................

Stress-stretch behavior and volume-stretch behavior for equibiaxial
loading of a saturated gel with M = 100 and X = 0.425 showing
the effect of different é. .........................

Stress-stretch behavior for equibiaxial loading with M = 100, X =
0.425, and 6:0. Three nonsaturated response curves are depicted.
Each nonsaturated response branches off of the saturated response
“backbone curve”. ............................

Stress-stretch behavior and volume-stretch behavior for uniaxial load-
ing of a saturated gel with M = 100 and X = 0.425 showing the effect
of different 5. ...............................

Stress-stretch behavior for uniaxial loading with M = 100, X = 0.425,
and .5 =0. Three nonsaturated response curves are depicted .......

The dependence of volume and stress on uniaxial stretch for M =
100, X = 0.425, and 6 = 1. Three nonsaturated examples are shown
corresponding to loss of saturation at )‘fmi— A = 15,16 and 1.7. . . .

Stress-stretch behavior and volume-stress behavior for uniaxial testing
of a saturated gel as given by R. Monroe in [74] .............

viii

20

25

32

33

36

38

40

42

3.7

3.8

3.9

3.10

4.1

4.3

4.4

4.6

Stress-stretch behavior for equitriaxial loading of a saturated gel with
M = 100 and X = 0.425 showing the effect of different 5. For g > 0.038
the graphs are monotonic. For 0 < 5 < 0.038 the graphs are decreasing

on a finite interval in g ...........................

Critical value Ecrit for the existence of a local stress maximum in the
equitriaxial stress response. For 6 > gem-t the response is monotonic. .

Stress-stretch behavior for equitriaxial loading with M = 100, X =
0.425, and {=05 showing three examples of nonsaturated response.
Since the stretch A = J1/3, no additional stretch is possible after loss
of saturation. ...............................

Comparison of a saturated and a nonsaturated gel under uniaxial load.
On the left is a comparison of the hard device energy expression. On
the right is the resultant force as a function of stretch. As in Figure
3.5 the constitutive parameters are given by M = 100, X = 0.425 and
5 = 1. The transitions are taken to occur at A = 15,16 and 1.7. . . .

Radial deformation (1‘ = 'r/ ((1%)) and local volume change ratio as
a function of radial position. As in previous figures, M = 100 and
X = 0.425. The cylinder geometry is such that R0 = 2R,- ........

Resultant axial force as a function of mechanical stretch and total
volume as a function of mechanical stretch for both the everted cylinder
and the uneverted cylinder using the same material parameters as in
Figure 4.1. For the everted cylinder, the reference geometry is again
such that R0 = 2B,. ...........................

Total volume after eversion as a function of mechanical stretch, and
resultant force as a function of mechanical stretch, for various values
of 5, again using M = 100, X = 0.425 and R0 = 212,-. .........

Local volume change J versus A; for an everted cylinder with R0 = 2R,-
and material parameters M = 100, X = 0.425 and either E = 0 or
6 = 1. On the left, .5 = 0 and the loss of saturation corresponds to
that depicted in Figure 4.5. After loss of saturation, increasing stretch
redistributes the finite amount of fluid from the outer region to the in-
ner region. On the right, E = 1 and the loss of saturation corresponds
to that depicted in Figure 4.7. For the nonsaturated cylinder, increas-
ing stretch now redistributes the finite amount of fluid from the inner
region to the outer region. ........................

44

45

50

63

69

74

4.8

5.1

5.2

5.3
5.6

5.8

5.9

5.10

5.11

Resultant axial force as a function of mechanical stretch for the everted
cylinder with R0 = 2B,. Material parameters correspond to those in
Figure 4.3(c), namely: M = 100, X = 0.425, and f = 1. Three different
values of qu are considered, corresponding to A: = 0.8, 0.9, and 1.0.
The nonsaturated response curve is first above, and then below, the
saturated response curve, and so mirrors the case of this material under
homogeneous deformation (Figure 3.10(b)). ..............

Free swelling volume change J f3 and stretch C for different a on the
basis of (2.3) and (5.1) using (5.3). ...................

J normalized by the free swelling value J f3 = C3 as a function of
R = R/R, (panels (a) and (c)); and the twist g = 6 — 6 as a function
of R (panels (b) and (d)), for boundary conditions given by (5.19) and
(5.20) with Ro/R1 = 2. The problem is equivalent to that presented
in [48] and panels (a), (b), (c) correspond respectively to Figures 1, 4,
5 in [48] ...................................

Normal stress 01-,- for the same problem as described in Figure 5.2.

Additional graphs for the boundary value problem described in Figure
5.3. The dashed lines in panels (a) and (0) give the baseline asso-
ciated with a pure polymer (no liquid) by showing the volume ratio

_, —-1
Vpoly/st .. J f3 ..............................

Additional graphs for the boundary value problem described in Figure
5.4 ......................................

0,»,- (Ti) versus 6 for a = 0. Panel (a) is the magnified portion of panel
(b) near 6 = 1. Panel (0) provides additional curves near 6 = 1 with
the top curve for T” = 0.8; bottom curve for T * = 1.0, and with
additional curves for T* increments of 0.01 in between. ........

(a) Torque vs. twist relation with a = 0 for the boundary conditions
(5.24). Here T* = 0.937, 21) = 0.90 corresponds to the local maximum.
Also T * = 0.377, 1,0 = 1.70 is the effective limit of our numerical
calculation and corresponds to 6 = 0.15. (b) Volume ratio V/st
versus twist angle 1/2 for the same problem. ..............

The post-maximum solutions for a = 0. These correspond to a con-
tinuation of panels (a) - (d) of Figure 5.4. Now T* decreases with
increasing 1,0 as is evident from panel (d) .................

78

82

90

91

97

100

103

105

106

5.12

5.13

5.14

5.15

5.16

6.1
6.2
6.4

6.5

6.6

Mechanical behavior for the nonsaturated solutions when loss of satu-
ration takes place at a standard solution with T* = 0.6. Here a = 0
and Vh-q = 40.469LRg’. The solid curve in each panel corresponds to
the instant of transition while the other curves correspond to increas-
ing T*. The values of 21) associated with T* = 0.60, 0.80, 1.00, 1.20,
1.40 are u ==0.418, 0.562, 0.708, 0.858, 1.011 respectively. .......

Mechanical behavior for nonsaturated solutions when loss of saturation
takes place at the post-maximum solution with T * = 0.6. Here a =--- 0
and Vliq = 53.371LRg. The solid curve in each panel corresponds to
the instant of transition while the other curves correspond to increasing
T*. The values of 1p associated with T* = 0.60, 0.80, 1.00, 1.20, 1.40
are «p = 1.412, 2.053, 2.800, 3.632, 4.533 respectively. .........

Saturated and nonsaturated torque-twist response for a = 0. Two
examples are shown, each corresponding to T * = 0.6. The square
markers on the nonsaturated curves represent computational points
and the intervening dots are to guide the eye. The transition on the left
occurs at a standard solution (corresponding to the solutions displayed
in Figure 5.12). The transition on the right occurs at a post-maximum
solution (corresponding to the solutions displayed in Figure 5.13). . .

Saturated, transition and nonsaturated states associated with loss of
saturation taking place from a standard solution when 1,0 = 0.42. Thus
if = 0.62 corresponds to a nonsaturated state. For comparison, the
response for a hypothetical ’t/J = 0.62 saturated state is also shown. . .

Saturated, transition and nonsaturated states associated with loss of
saturation taking place from a post-maximum solution when 2,!) = 1.41.
Thus 2,!) = 1.60 corresponds to a nonsaturated state. For comparison,

the response for a hypothetical w = 1.60 saturated state is also shown.

Mapping of fluid and solid particles ...................
Free swelling of a fiber reinforced layer .................

The variation of the change in orientation angle w — (2 (in radians)
with g = o and M/p = 100, x = 0.425. .................

Representation of a pressure differential that will cause diffusional fluid
seepage through a fiber reinforced layer .................

Comparison of the numerical results with the experimental data [24] .

110

112

113

114

116

120

127

135

136

142

6.7

6.8

6.9

6.10

6.11

6.12

6.13

6.14

6.15

6.16

6.17

The variation of nondimensional displacement in the thickness direc-
tion after free swelling at PH — P0 = 200psz' with 'y = 0 for C = 0 and

5:1. ...................................

The normal displacement of the free surface vs. pressure difference for
the material parameters in Figure 6.7(a) with C = 0 and Figure 6.7(b)
with C = 1 ................................

Variation of the fluid flux q/pg and the required reactive shear stress
07,3,”(11) with pressure difference PH — P0 for various fiber reinforce-
ment coefficients ‘7 with C = 0, n = 5, a = 1.556 x 1019 gm/cc - day

when the boundary conditions are given by (6.67) ............

The fibers become curved in the inhomogeneous deformation associated
with pressure driven fluid seepage as depicted in panel (c). This is a
schematic representation ..........................

The fiber geometry at various seepage rates with 7y = 2.0, C = 0, n = 5
and a = 1.556 x 1019 gm/ cc - day when the boundary conditions are
given by (6.67). (xf denotes the projected distance onto x-axis between
an arbitrary point on the fiber and the tip (z=0) of the same fiber.) .

Variation of the fluid flux q/pg)c and shear displacement f (H) — I‘(H)
with pressure difference PH - P0 for various fiber reinforcement coef-
ficients 7 with C = 0, n = 5, a = 1.556 x 1019 gm/cc - day when the
boundary conditions are given by (6.68). ................

The shear displacements as a function of the initial fiber orientation
under a free surface traction boundary condition at Z = H .......

The normal displacements as a function of initial fiber orientations
with a free surface traction boundary condition at Z = H. ......

The fiber geometry at various seepage rates with 7y = 2.0, C = 0, n = 5
and a = 1.556 x 1019 gm/cc - day when the boundary conditions are
given by (6.68). The common slope value at Z = H in (b) is a direct
consequence of the boundary conditions (6.68)4 and (6.68) 5. .....

An explicit forward time-stepping scheme ................

Change in compaction during transient diffusion in the absence of fibers
('7 = 0) for a material with C = 0, M / p = 100. In this simulation the
steady state seepage with PH —— P0 = 501932" is abruptly altered to
PH -- P0 = ZOOpS’i. ...........................

xii

143

144

146

147

148

149

151

152

153

157

6.18 The variation of pf/pg along the thickness for the simulation consid-
ered in Figure 6.17 .............................

6.19 The variation of normal stresses along the thickness direction for the
simulation considered in Figure 6.17. ..................

6.20 The amount of fluid in the layer changes with time: from its initial
steady state value (where PH — P0 = 50psi), and finally approaching
the final steady state value (where PH — P0 = 200psi). .......

6.21 The change of interstitial fluid amount with time for various pressure
difference. Free swelling state (PH — P0 = 0) is abruptly altered when
PH — P0 is suddenly changed to 50psi, 100psi, 200psi, 300psi, and
400psi, respectively .............................

6.22 The shear and normal contraction associated with fiber reinforced tran-
sient diffusion where '7 = 1.0, Q = 45°, and C = 0. ...........

6.23 The shear and normal contraction associated with fiber reinforced tran-
sient diffusion where ’7 = 10.0, $2 = 45°, and C = 0. ..........

xiii

160

161

163

166

Chapter 1

Introduction

By a gel we mean a system of crosslinked polymer chains mixed together with
a low molecular weight liquid. In the liquid environment, gels could undergo a large
deformation associated with swelling or shrinking. Their wide applications in drug
delivery [1,2], surgical dressings [3], artificial tissue [4] (skin, articular cartilage, etc.)
and control material (micro-control [5], smart materials, etc.) in engineering motivate
the desire to better understand their underlying mechanical response. Depending on
the nature of the elastomer and fluid, a variety of ambient conditions, such as tem-
perature [6,7], pH of the liquid bath [8,9], and light exposure [10, 11], could possibly
determine the amount of liquid that is absorbed into or expelled out of the gel. This
is typically coupled with the mechanical deformation of the polymer network. Ther-
modynamical and mechanical balance principles describe these processes, in which
case a central issue is the development of a modeling framework that is sufficiently
general to accommodate these effects upon specification of appropriate constitutive
relations. One possible framework is that of small deformation; such theories describe
swelling by allowing for parametric dependence of constitutive parameters (e. g. elastic
constants, osmotic pressure) on the underlying electrochemistry (e.g. concentrations,

charge density) [12]. Such linear descriptions are useful for describing mechanical

behavior in fixed ranges wherein issues of large strain, instability, and transition from
saturation to unsaturation do not arise. However, in certain situations such phenom-
ena cannot be avoided, or will even be exploited, in which case it is necessary to

consider a finite strain description.

1.1 Background

Gels are high molecular weight crosslinked polymer networks immersed in a low
molecular weight liquid. When water is the liquid, then these are typically referred
to as hydrogels. The liquid can also be organic in nature as for example discussed
in Chapter 7 of TIeloar’s well known treatise The Physics of Rubber Elasticity [13].
Although hydrogels would generally be much softer than a rubbery material infused
with an organic liquid, both are regarded as gels for the purpose of this thesis, and
we shall use the term gel throughout in this broad sense.

Some of the earliest experiments on the stress-stretch behavior of swollen gels
were performed by Gee [14]. Simple uniaxial tension tests were made on a range of
vulcanized rubbers at various degrees of swelling. The experimental results depart
greatly from those predicted by the standard large strain theories of rubber elasticity.
These classical theories are associated with a Gaussian statistical mechanics treatment
of the polymer chains. Similar tests can also be found in [15] and [16].

McKenna et al. [17] carried out torsional and uniaxial compression tests on both
dry and swollen cylindrical rubbers. They concluded that the dry-state elastic free
energy function could also describe the elastic behavior of swollen gels for non-highly
cross—linked rubbers. This implies that the presence of liquid molecules does not affect
the elastic component of the free energy function. Their experiments are claimed to be
consistent with the Frenkel-Flory-Rehner assumption that the total energy of swollen

networks is comprised of the free energy of mixing and the elastic free energy.

There are also many models [12,18] that treat equilibrium mechanical behavior
of gels in the setting of the classical small strain tensor 8 = %(Vu + VuT) (where u
is the displacement vector) within the theory of linear elasticity. These include the
small strain linear theory of poroelasticity [19]. They are not useful for our purposes
and will be omitted from what follows. Instead, it is necessary to use large strain
elasticity theory.

A separate, and complicating issue is the internal flow of fluid in gels so as to
achieve equilibrium if the system is not initially in equilibrium. The diffusion of fluid
through polymer networks has been investigated a great deal [20—23]. One of the
most frequently referred to modern set of experiments on this can be found in [24],
where Paul and Ebra-Lima measured the non-linear diffusion of twelve organic liquids
through a swollen gel membrane. The diffusion was induced by the pressure difference
across the membrane. The diffusion coefficients and also the material parameters
of both the fluid and the hyperelastic rubber material were determined from their
experimental results. These material constants were then widely adopted in later
modeling frameworks.

One of the earliest continuum theories that can be used to describe the large
strain behavior of incompressible hyperelastic materials was presented by Mooney [25]
and Rivlin [26]. This apples to rubber-like materials but does not, in its original
presentation, account for gel mixtures where the interaction of the elastomer and
liquid is present. Their basic theory involves an energy argument that the elastic free
energy (D depends only on the three principal stretches A1, A2 and A3. In current

notation, their model is equivalent to the expression

«p = <1>(A1,A2,A3) = g [(1 — g) (A'f‘ + A3 + A3 — 3) +6 (£93 + AgAg + A§A§ - 3)]
(1.1)

where p > 0 is the shear modulus for infinitesimal strains and C (obeying 0 g C S 1) is

u

the Mooney-ijlin adjustable parameter. Because this model is for an incompressible
hyperelastic material, implying volume preservation during mechanical deformation,
it leads to a constraint condition A1A2A3 = 1 on the three stretches. In the case

associated with C = 0, eq. (1.1) retrieves the neo-Hookean strain energy function.
,u
<I>=<I>(A1,A2,A3)= E (A§+A§+A§—3), (1.2)

which arises from a Gaussian statistical treatment of polymer chain entropy.

One of the basic issues for the mechanical behavior of gels is how much they
swell depending on their liquid environment. In the absence of load this is called free
swelling. For a given value of free swelling, the question arises as to how an addi-
tional applied load on the swollen gel further deforms the gel. Various hyperelastic
models have been used for this purpose. In other words, these models seek to predict
how stresses applied to a swollen gel cause additional stretch beyond that associated
with free swelling. The more chemistry oriented literature often accomplished this
by balancing the elastic free energy (as a function of the stretches) with the work
done by the principal stresses 0‘1, 02 and 03. This then gives the relation between
the principal stresses and the principal stretches. In addition to the Mooney-ijlin
material mentioned above, some of the other energy functions that have been used for
this purpose include the constrained chain model [27, 28], localization model [29,30],
liquid-like model [31,32], and eight-chain model [33]. Han et al. [34] compared these
models and concluded that the constrained chain model fits experimental data [35—38]
best for both dry and swollen states, whereas the eight-chain model of Arruda and
Boyce shows less agreement with the experimental data. However, later Arruda and
Boyce [39] questioned the calculation executed by Han et al. [34] and put forward
their eight-chain non-Gaussian Model in order to account for the same trends as

in the experimental data for large stretches. On this basis, they then proposed a

hybrid model that suitably combines the F lory-Erman constrained chain model (so
as to capture the main features of stress-strain behavior for small strains) and the
Arruda-Boyce eight-chain model (which was developed to describe the large strain
behavior well). The strain energy density function of this hybrid model gives rise to

the following stress-strain relationship for uniaxial loading:

1
2 ‘3
Nke n _ V Ac 1 1 I
VET-£1: 1 375— (A2 — A) + ENJkez/p (A1A2 — Azx)

(1.3)
Here up is the volume fraction of the polymer relative to the dry state and represents
the degree of swelling, N is the number of chains in the polymer per unit unswollen
volume, N J is the junction density, lc is Boltzmann’s constant, 9 is absolute temper-
ature, n is the number of links in the chain, Ac = [(A2 + 2/A)/3]1/3, £—1( ) is the
inverse Langevin function, and A1, A2 are both functions of A and up. These are too
complicated to be discussed in detail here and the reader is referred to [39] for more
information.

In summary, the treatments described above in [27—34] seek to use the theory of
hyperelasticity to model gels after free swelling has occurred. They account for the
strain energy function of hyperelastic materials and give a reasonable linear approxi-
mation to the neo-Hookean model. The main distinction between different models lies
in the modification of the strain energy function for moderate and large deformation.

In order for a model to also determine the amount of swelling, it is necessary to
connect the mechanical response to the interaction between the polymer network and
the liquid which together make up the gel. The attempt to make this connection was
made over half a century ago. The swelling of the gels depends upon the interaction
between the solid elastomer and the surrounding liquid. One of the earliest exper-
imental measurements on this interaction goes back to [40], where Gee and Treolar

gave a thermodynamic description of the rubber-benzene mixture system by measur-

T

ing the vapor pressures. A statistical mechanical theory of the equilibrium degree of
swelling was presented by Flory [41] and Huggins [42,43] separately. Based on this

theory, the free energy of mixing H is
H = M[(u-1 — 1)ln(1 — u) + X(1— 11)], (1.4)

where M and X are material parameters, and V is the volume fraction of the elastomer
in the mixture. They seek to determine the free swelling ratio 1/1/ by an energy
argument that seeks to account for: (1) the mixing entropy of elastomer and liquid; (2)
mixing enthalpy of elastomer and liquid; and (3) strain energy of the elastomer. If a
system is originally in non-equilibrium, then it is assumed that internal fluid diffusion
will occur until equilibrium is achieved. This diffusion process can be described by a
mixture theory approach. For now we restrict attention to equilibrium configurations
once any fluid redistribution is complete.

On the basis of Flory-Huggins theory, Tfeloar [44] more formally introduced an
extra term that explicitly accounts the work due to tractions when load is applied.
By this treatment, the swelling equilibrium associated with homogeneous deformation
under simple loading conditions could be determined. Treloar [45] later considered
the problem of a cylinder subjected to combined axial extension and torsion about
the axis. For large deformation the radial distribution of the stress, strain, and the
amount of swelling are obtained by numerical analysis. It is found that the gel volume
decreases with the torsion.

Wineman and Rajagopal and co—workers constructed a continuum mechanical
treatment incorporating finite strain on the basis of mixture theory, as discussed in a
series of papers [46—48]. In [48] they consider a static problem where a relative rotation
was applied to the two lateral surfaces of a swollen hollow cylinder. Their main

purpose was to find the fluid redistribution with various twist moments. However,

1'

they limited attention to boundary conditions that did not allow for volume change
after swelling. In an even earlier paper [49], Wineman, Rajagopal, and co-workers
consider a torsion problem that does treat load dependent volume change using a full
mixture theory formulation. As in [45], they find that torsion causes volume decrease;
this corresponds to the expelling of fluid.

In this work, I have examined similar static boundary value problems with dif-
ferent loading conditions. Both volume decrease and volume increase associated with
the applied load can be obtained. The latter case corresponds to fluid absorption.
For these static problems it is not necessary to involve mixture theory and so are
similar to the classical approach of Tmloar [45]. For the case of fluid absorption, it
may further be the case that insufficient fluid is present to attain a saturated state.
In this case, the associated transition from saturation to unsaturation is considered.
In particular, it will be indicated how the overall mechanical response becomes stiffer
if this transition takes place will be indicated.

These are not the only continuum models that have been developed to treat the
swelling behavior of gels. Instead of introducing the mixing free energy, Marra and
Ramesh and collaborators [50,51] introduce an evolution internal variable a, which is
defined as the volume ratio of the interstitial fluid in the current state to that in the
reference state. It is assumed that this internal variable acts as the only independent
variable of the specified free energy function. To describe the volume change (called
actuation in their article) they consider the total free energy function as the sum of
elastic energy and an additional energy term due to the effects of free actuation, which
means uniform expansion or contraction of the polymer in the fluid. This corresponds

to the statement, in current notion

W : ©(11712) + g(I3,Ct), (15)

 

where W represents the total free energy and <I> can be regarded as some elastic energy
that is appropriate for an incompressible material. This (I) is typically taken to be
in the neo-Hookean, Mooney-Rivlin or even the more general Ogden [52] form. The
function g in (1.5) is claimed to account for the effects of free actuation. The three
invariants and the material parameters in <I> and g are all assumed to be functions
of a. Then an evolution law for the internal variable a is developed. Uniaxial and
equibiaxial loading tests were performed. Experimental data are obtained and used
in fitting the material parameters in Ogden energy function. These parameters are
then substituted into a finite element model. It was then claimed that the numerical
results show good agreement with the data.

Certain hydrogels exhibit phase changes of various kinds such as an abrupt vol-
ume increase when either the temperature or chemistry of the liquid is changed. These
are sometimes called stimulus-repsonse hydrogels. Dolbow et al. [53—55] presented a
continuum model for volume transitions in hydrogels that exhibit a sudden collapse in
volume at a specific temperature. Their model includes a sharp-interface separating
swollen and collapsed phases. With obeying the deformational and diffusion poten-
tial coherence across the interface, a set of equations of force balance and chemical
potential balance were derived in the domain of each phase and on the interface as
well.

Another aspect of interests in the swelling behavior of gels is the study of fluid
diffusion through hyperelastic material. A variety of theoretical treatments [56, 57]
have been developed to model the mechanical and thermodynamical behaviors of such
diffusion process. A common-used biphasic theory called the theory of porus media
(TPM), involving the mixture theory and volume fractions, is applicable to such non—
linear diffusion problems. One of the earliest presentations of this treatment can
be found in the work of Bowen [58]. His framework, which involves partial stresses

for the solid and fluid phases, makes use of mass balances, momentum balances, an

energy balance law and the Clausius-Duhem inequality.

Many material models generalizations were then subsequently embedded into this
binary mixture theory, such as a viscoelastic model [59] and an elasto-viscoplastic
model [60, 61]. Different numerical methods based on finite element analysis have
also been developed and used to study the phenomenon of consolidation in the porous
media, including Galerkin finite element method and least—square mixed finite element
method [62], [63]. Within this framework, Markert [64] has presented a 3-D finite
element analysis on the nonlinear fluid flow through a porous polyurethane form.

Ehlers et al. [65,66] followed a related energy and entropy balance treatment
and used it to simulate the elastic deformation of liquid saturated porous solids. In
particular these works introduce the idea of a fully compact material meaning that
all the internal pores have been closed, squeezing out all the fluid. Their energy
functions are designed to give this full compaction only in the limit as appropriate
normal stress components trend to negative infinity.

The compressibility condition on the solid and the fluid phase of the mixture
is always an important fact to be considered in TPM. Palomar and Doblare [67]
utilized an Augmented Lagrangian formulation to enforce the incompressibility con-
dition on both the solid and fluid phase based on a fibre-reinforced porohyperelastic
model. Diebels extended the mixture theory and applied it to both incompressible
and compressible binary systems [68].

Specifically, [61, 67, 69] deal with material anisotropy. The numerical examples
in [67] show how the additional stiffness due to the fiber reinforcing causes the fibrous
materials to exhibit less deformation under the same amount of extensile load. Among
the more interesting results obtained in [67] is the demonstration that the introduction
of fiber network leads to a higher relative fluid diffusion velocity as well as a higher
pore pressure in an unconfined compression simulation.

Under the TPM framework, Markert et al. [70] derived a set of full contact

 

boundary conditions at the interface of two saturated porous media with the assump-
tion that there is no dissipative effects at the interface. These conditions include the
solid velocity continuity condition, the normal seepage velocity continuity condition,
the pore-fluid pressure jump condition and the solid effective normal stress jump
condition. The latter two conditions have the property that if two identical porous
bodies are put in contact with each other then both the pore pressure and the solid
effective normal strass are continuous. However if the two contacting porous bodies
have different solid-fluid phase ratio at the interface, then the pore pressure and the
solid effective normal stress exhibit a jump and mixture theory arguments are used
to obtain detailed expressions for these jump (see equation (47) and (48) of [70]).

However none of the works [56—70] obviously takes into account the interacting
behavior between the fluid and the solid skeleton. To also connect the mechanical
response to the interaction between the polymer network and the interstitial fluid,
Rajagopal, Wineman and collaborators (Rajagopal et al. [46], Gandhi et al. [47] and
references therein) introduced a Flory-Huggins type term in the free energy expression
so as to account for the interacting continua in the diffusion problems. Pressure-
induced diffusion through spherical and cylindrical geometries were studied in [46,47].
Diffusion through a porous swollen layer with lateral stretch and shear deformation
was also considered in [47]. In particular these works provide numerical solutions to
boundary value problems with one lateral surface held fixed while only the pressure
was specified on the other surface.

These treatments also pointed out some controversies with respect to specify-
ing the tractions on the solid and fluid boundaries. In an alternative but related
variational treatment, Back and Srinivasa [71] extremize the Helmholtz potential of
the system containing the swollen solid and the fluid over all admissible system con-
figurations. This variational procedure involves the mass conservation equations for

both the solid / fluid mixture phase and the surrounding pure fluid phase. This results

10

in a set of governing equations and boundary conditions. Energy dissipation is also
considered in [71] to account for frictional effects as the fluid diffuses through the
elastomer.

Other means for arriving at this framework are also possible. This includes that
given recently in [72] and [73]. In partcular, the treatment in [72] introduces the the
liquid phase concentration as an independent variable prior to the requirement that
the liquid component fraction in the gel be given by 1 — V. The equivalent of this
constraint is enforced by a Lagrange multiplier in [72]. The condition for free swelling
is then found to be given by a requirement that the chemical potential for the pure
liquid, say ”liq: is equal to an appropriate variation of the gel’s total free energy with

respect to a change in concentration.

1.2 Objectives and Scope

The main purpose of this research is to model the mechanical deformation and
also the redistribution of the concentration of fluid throughout the system in response
to changes in load or displacement at the system boundary. In so doing, it distin-
guishes between saturated and nonsaturated gels. Numerical analysis is introduced
to obtain the equilibrium deformation of the polymer network, from which the redis-
tribution of the concentration of fluid is derived. This is based on the assumption of
volume additivity that the volume of the mixture is the sum of that of the elastomer
(solid phase) and that of the interistic fluid (fluid phase).

In particular, even though both the polymer matrix component and the fluid
component of the gel system are regarded as individually incompressible when iso—
lated from the other component, the overall gel system, as long as it is saturated,
is then described in terms of a hyperelastic framework in which the gel is nomi-

nally compressible. This is because a fixed amount of polymer component can be

11

T

regarded as defining a gel patch (or gel “element”). The amount of fluid component
within such a gel element can change, giving rise to a change in the volume of the
gel element. In other words, redistribution of the fluid component with respect to
the matrix component generates a local volume change. This is in keeping with the
Treloar’s treatment. He writes in [13] that if a rubber that is incompressible in the

absence of a liquid swelling agent is subsequently swollen, then

“ the swollen rubber in continuous equilibrium with a surrounding
liquid may be regarded, from the purely formal standpoint, as having

mechanical properties equivalent to those of a compressible material. ”

In particular, since the rubber that is referred to in the above quote is regarded as
being in contact with a liquid bath, the situation refers to a state of liquid saturation.
Thus the saturated equilibrium response is described by a single stored energy density
function W. The present work has taken the Frenkel-Flory-Rehner assumption that
total energy of the gel is comprised of elastic free energy and mixing free energy, which
are additive. The choice of the specific form of elastic free energy is not the central
problem in this work, therefore the specific Mooney-Rivlin model will be adopted
for modeling isotropic gels. The Flory-Huggins equation (1.4) will be used for the
mixing free energy. Then the total free energy W for a swollen gel depends only on
the deformation gradient F through the principal stretches A1, A2, A3, or equivalently
the three invariants 11, 12, I3 of the left Cauchy-Green deformation tensor B, namely

FFT. By virtue of (1.1) and (1.4), this becomes
W = <I>(A1, A2, A3) + H(V). (1.6)

The saturated Cauchy stress tensor then follows from this stored energy density as

12

 

in compressible hyperelasticity
a = ——FT. (1.7)

The local volume change in the deformation then directly correlates to the local
change in fluid content.

If, however, the amount of liquid is limited, and in particular assume that all
the liquid is imbibed at some point in the loading process, the system then becomes
nonsaturated. A constant reaction stress —pI is generated so as to enforce the global
constraint that the overall volume of the gel remains fixed at the instant when the
transition from a saturated state to a nonsaturated state occurs. This gives the

nonsaturated Cauchy stress tensor
= ————FT — pI. (1.8)

Thus (1.7) and (1.8) give the relation between stress and deformation for homo-
geneous deformation (F independent of location). For nonhomogeneous deformation

without any body forces the stress field must satisfy
div 0 = 0. (1.9)

These equations are usually nonlinear and hard to obtain analytical solutions. There-
fore it is necessary to invoke numerical routines in order to obtain simulated results.

The above hyperelastic framework could describe equilibrium situations in which
both the fluid component and the polymer matrix component of the system are at
rest. In particular, each material point in the gel is regarded as a two component
mixture of polymer matrix and interpenetrating liquid. For a sort of more complicate

problems where there is relative motion between the fluid and matrix components, so

13

it is necessary to invoke the broader continuum mechanical framework that specifically
deals with separate mechanical balance principles for each component. This broader
framework is known by a number of names, including: large deformation mixture
theory, the theory of interacting continua, and the large deformation biphasic theory.
To this end, the process of fluid diffusion through hyperelastic media is investigated
in this work, where a major focus is on time dependent swelling as liquid diffuses
within the elastomeric matrix.

The structure of this work is outlined as follows. The hyperelastic theory for the
static mechanical response of swollen gels is developed in Chapter 2. The discussion
of the stress-strain swelling response under homogeneous deformation is discussed in
Chapter 3. Boundary value problem of an everted tube subject to an axial load is
considered in Chapter 4. The inhomogeneous deformation response of the axially
loaded tube is compared to the homogeneous deformation response of the axially
loaded tube when it is not everted. Chapter 5 deals with a boundary value problem of
a hollow tube with a relative twist at its lateral surfaces. Different radial displacement
and traction boundary conditions are considered along with the prescribed twist. In
Chapter 6, certain flow problems where fluid exhibits both steady-state and non-
steady-state seepage through the hyperelastic gel are studied. These seepage problems
consider the possibility of pressure driven fluid diffusion through a fiber reinforced
gel, such that the isotropic diffusion problem is retrieved as a special case.

For the equilibrium mechanical response discussed in Chapter 2 through Chapter
5, it is a main object to investigate the effects of fluid saturation and loss of saturation
in the system on the polymer network deformation as well as fluid redistribution
under certain mechanical loadings with different boundary conditions. While as for
fluid seepage problems it is assumed that there is always enough fluid flowing through
the polymer network so that there is no necessity to consider the loss of saturation

at all.

14

Chapter 2

Modeling of the Swelling Behavior
of Elastomeric Gels

2.1 The Flory-Huggins Equation for the Determi-
nation of Free Swelling

The most basic description of large strain effects in elastomeric gels concerns the
determination of the amount of fluid that perfuses a polymer matrix when it is placed
in a liquid bath. Let V denote the volume fraction of polymer matrix so that 1—1/ is the
volume fraction of the fluid component. A standard development proceeds by taking
molecular chain arguments for configurational entropy of crosslinked macromolecules
within a liquid bath, and coupling these to a phenomenological description of enthalpy

of mixing [13,44]. The requirement of a stationary free energy then leads to the

 

equation
M 1 1 2 1/3 — 0 21
n ’st +st+Xst + pus — (.)
x v 1 \_ _,
dilution crosslinking

for the free-swelling value of V which, as indicated above, will be denoted by ”far Here
M > 0, x and p > O are constitutive parameters that may vary with temperature,
pH, fluid chemistry, and other environmental factors. As indicated in (2.1), two
separate effects can be identified, that associated with dilution and that associated

with crosslinking. In the absence of crosslinking, the equation reduces to that for

15

the effect of dilution alone, and this is the equation that is usually referred to as
the Flory-Huggins equation [13]. The dilution term itself models the two effects of
mixing entropy and mixing enthalpy. The term containing x models the latter and is
positive for the standard case wherein polymer-polymer and liquid—liquid grouping is
favored over polymer-liquid grouping. In particular, larger x favors polymer-polymer
aggregation so that st increases with x. The remaining part of the dilution term
models the entropy contribution of polymer-liquid mixing to the free energy. In
particular, the parameter M is identified as the product of the ideal gas constant and
the absolute temperature divided by the molar volume of the liquid.

The utility of a theory that delivers (2.1) is due in no small part to the fact that
(2.1) has a unique solution st that takes values on the physically relevant interval
0 < st < 1. This solution is here referred to as the free swelling polymer volume
fmction. It is completely determined by x and the ratio u/M. Thus one may then
write ”f3 = 19f3(x, u/M). Increasing either X or u/M favors a more tightly bound
gel and so increases 19f3(x,p/M).

The connection between free swelling and hyperelasticity follows by introduc-
ing the deformation mapping that takes the polymer component from its original
unswollen location X to its stressed, swollen location y. Let J = detF where
F = 6y/6X is the deformation gradient of the mapping y(X). Correspondence
between J and u is immediate if both the polymer component and the fluid compo-
nent are individually incompressible since simple mixing then gives the identification
J = 1 / V. Thus J is the volumetric swelling of the gel as measured with respect to

the polymer component before liquid was present. It is therefore required that
J 2 1. (2.2)

Historically, the original development leading to (2.1) employed a notion of physical

16

—t_

variations. In particular, rather than using physical arguments to obtain an explicit
energy to minimize, the original papers of Flory and 'Ireloar employ an argument
which requires that there is no change in the free energy of the system in the event
that a single fluid particle is transferred between the gel and the pure liquid which
surrounds the gel. Although the free swelling volume fraction appears in (2.1), the
original arguments leading to this equation actually depend more on the concept of
occupied volume rather than the related concept of volume fraction. The connection
between occupied volume and volume fraction then follows from J = 1 / V. The
argument involving fluid particle transfer can thus be viewed as a variation with
respect to J. Accordingly, the energy that is formally minimized can be found by
substituting st —> 1 / J in (2.1), integrating the resulting expression with respect to
J, and, if desired, returning to the original volume fraction variable via J —+ 1/uf3.
This gives, to within an integration constant, an energy expression with the following

terms

1 3p —2/3
—MXst +M(V—f;-1)ID(1—st)+ —2—'st

effects of mixing entropy

\

effeCtS 0f mixing enthalpy effects of network elasticity
(2.3)

Hence (2.1) is equivalent to determining the free-swelling value of J as the root of

gag-(HM) + mm) = 0, (2.4)

with
H(J) = M((J — 1)ln(1 —- J—1)+ x(1- J—1)), (2.5)

and
\IJ(J) = 37“.!2/ 3. (2.6)

Here H + \II is the overall free energy, the former accounting for polymer-fluid inter-

17

action (mixing entropy and mixing enthalpy), the latter accounting for elastomeric
deformation of the crosslinking network. The limiting case a —+ 0 corresponds to the
absence of crosslinking and hence no elastic effect. Equation (2.4) then reduces to

h(J) = 0 where we have introduced the notation
h(J) E —. (2.7)

It is to be remarked that (2.2) and (2.4) can be regarded as a general framework, which
retrieves the well known equation (2.1) once the specific mathematical forms (2.5) and
(2.6) are invoked. Thus statements involving the functions H and \Il need not be tied
to (2.5) and (2.6) unless specifically indicated. Similarly, (2.7) is a definition for h
that is not tied to any constitutive form. For the specific constitutive function (2.5),
h is given by

h(J) = M[1n (1 — J‘l) + J‘1 + XJ—2]. (2.8)

The basic behavior of this framework with the usual forms (2.5) and (2.6) is sensitive
to whether X < 1/2 or X > 1/2 as described next.

If (2.5) and (2.6) hold and X > 0.5 then there is a unique value of J > 1
that causes h to vanish. This value of J corresponds to free swelling in the ab-
sence of crosslinking and we denote this value by Jnd where the subscript refers to
nogrosslinking. Thus h(Jncl) = 0 and Jncl = 1/19fs(X,0). Moreover h(J) < 0 for
1 < J < Jnd and h(J) > 0 for J > Jncl- Consequently Jnd provides a local minima
of H. One also finds that Jnd < 2X/(2X — 1) and that Jncl —> 00 as X —> 0.5.

If (2.5) and (2.6) hold and X < 0.5 then there are no solutions to h(J) = 0
because h(J) < 0 for all J > 1. Hence for X < 0.5 the infimum of H is achieved as
J —+ 00. Since h(J) —> O as J —» 00 it is convenient to formally define Jnd E 00

for X < 0.5. The derivative h' (J) > O for all J > 1. These qualitative features are

18

 

 

 

 

 

 

 

 

    
 

 

 

 

J
(a) X > 0.5
0.001 .....................
0.000
30.001
.1: \
— 0.002 \ X increases from 0.1 .
to 0.2, 0.3, 0.4, 0.45.
-0003 .......
20 40 60 80 100
J
(b) X < 0.5

Figure 2.1: The constitutive function h(J) in (2.8) for various x. Vaues of J giving
h(J) = 0 model gel dilation on the basis of mixing entropy and mixing enthalpy, but
do not account for the elastic effect of polymer crosslinking.

summarized in Figure 2.1.

Accounting for elastic interconnection (a > 0), equations (2.4) — (2.6) become
h(J) + uJT1/3 = 0. In conjunction with (2.8) this recovers (2.1) for the free swelling
value st E l/st. There is a unique finite st solution to (2.4) ~ (2.6) for all X
whenever p > 0, which is formally given by st = l/fif3(X,u/M). This solution
obeys the inequality J f3 < Jncl as would be expected since the crosslinking inhibits
swelling. Euther h(J ) +pJ _1/ 3 is respectively positive or negative as J is respectively

greater than or lass than J fs' Figure 2.2 indicates these relations.

 

T r r v *1

4.

     

h(J)/p+J‘1/3
O

\
X increases from 0.1 to 0.25, 0.45,
0.55, 0.75, 0.9

 

l L n | J l

 

 

0' 5 10111‘15‘ 20
J

Figure 2.2: The constitutive expression h(J) + ,uJ-l/3 for various X when u =

0.01M . The new term uJ‘1/3 describes the effect of crosslinks. Values of J giving
h(J) + uJT1/3 = 0 describe free swelling.

A simple interpretation concerns polymer that has some original volume, say
Vpoly! when free of liquid. Placing this nominally dry polymer into a liquid bath
causes it to swell by uptake of fluid. The resulting gel occupiae new volume J Vpoly
where J minimizes the free energy H (J) + \II(J). For the energy forms (2.5) and
(2.6) this free energy involves the three previously mentioned effects, one of which

favors swelling (the configurational entropy), one of which opposes it (crosslinking),

20

and one of which could go either way (mixing enthalpy) although it opposes swelling
in the standard case X > 0. In all cases, crosslinking (u > 0) ensures that the swollen
volume remains finite. If crosslinks are not present (,u = 0) then the mixing enthalpy
must sufficiently favor same phase agglomeration (X > 0.5) for the swollen volume to
remain finite. However if crosslinking is not present and the mixture enthalpy is not
sufficiently conducive to this same phase agglomeration (p. = 0 and X < 0.5) then
J —+ 00 and the polymer goes into solution, dispersing itself throughout the fluid
bath.

Implicit in the above discussion is a requirement that sufficient liquid is available

for saturation of polymer with the liquid. This requires that
JVpoly S Vpoly + Vliq’ (2.9)

where Vliq is the original fluid volume prior to introduction of polymer. If (2.9)
holds then the gel swells to its energetically favored saturation value. However if J
as determined on the basis of (2.4) is larger than that permitted by (2.9) then the

swelling is limited by the availability of fluid to the value
J* E 1 + Vliq/Vpoly (2.10)

which is less than the free swelling value J fs determined on the basis of (2.4). The
gel is then no longer saturated. In what follows the use of a subscript * or a su-
perscript * will denote a value that demarcates a transition between saturation and

nonsaturation.

21

2.2 Hyperelastic Constitutive Theory

In free swelling, the deformation gradient of the mapping y(X) is F = J}é3I
where I is the identity tensor. In order to treat more general deformations due to
the effect of mechanical loading, we now consider a hyperelastic framework. Let

B = FFT and C = FTF and let 11, 12 and I3 be the associated principal scalar

invariants

11 = TraceB, 12 = 12 — Trace B2 , I = detB = J2.
1 3

[\DIH

The overall stored energy density will be expressed as the sum of elastic free energy

of the polymer network CI) and a mixing free energy H
W(F) = <1>(11,12,J) + H(J). (2.11)

The mixing energy H is the same as that discussed in the previous section and its
derivative will continue to be denoted by h (viz. equation (2.7)). The specific Flory-

Huggins form (2.5) for H will be used in the examples that follow.

2.2. 1 Free Swelling

In the absence of mechanical loading, the associated free swelling is determined

by minimizing W in the class of simple volumetric expansions F = J 1/ 31 so that

I 1 = 3J2/ 3 and 12 = 3J4/ 3. Thus J is determined on the basis of
d 2/3 4/3 _
dJ (cp(3J ,3J ,J) + H(J)) _ o. (2.12)

Comparison with (2.4) indicates that correspondence with the free swelling frame—

22

work described in the previous section will follow provided that
<I>(3J2/ 3, 3J4/ 3, J) = \1/(J) + constant.
Correspondence with the particular Flory-Huggins form (2.6) will obtain if
<I>(3J2/ 3, 3J4/ 3, J) = §2EJ2/ 3 + constant.

In what follows we shall use the Mooney-Rivlin form for the elastic energy (I) of the

polymer network

¢=¢(11712)=§[(1—€)(11‘3)+5(12-3)[a (035$ 1). (2.13)

The case 5 = 0 gives the neo—Hookean specialization. For any E obeying 0 S 5 S 1 in
(2.13), one finds that

<I>(3J2/3, 3J4/3, J) = -§2’i((1—5)J2/3 + gJ4/3 — 1)

so that the special neo—Hookean case of 5 = 0 retrieves the original equation (2.1)
upon taking J = J fs = 1/1/f3. For the more general Mooney-Rivlin form with 6 > 0

equation (2.1) is augmented so as to contain an additional term and so becomes
1 3 —1 3
M[ln (1 — st) + st + Xugs] + u[(1 — {)l/fé + 251/ 8/ ] = 0. (2.14)

Introduce the free swelling stretch ratio

(2 J1/3 ____ V—1/3

3 f3 , (2.15)

where st is the root of (2.14). This root now depends upon g in addition to X and

23

p/M. As such, the free swelling behavior under (2.14) is not as simply described as .
was the case for (2.1). Some indication of the effect of the new constitutive parameter
5 upon st can be obtained by taking representative values for X and u/M. To this
end we follow [48] where the equivalent of a neo—Hookean (I) is considered for the
modeling of a vulcanized rubber in contact with toluene as considered by Paul and
Ebra—Lima in [24]. The following values are taken in [48]: M = 2.379 X 108dyne/cm2,
X = 0.425 and u = 2.375 x106dyne/cm2. Note for these values that the dimensionless
parameter M = M / p = 100.17 x 100. The value of the free swelling stretch ratio C
as a function of 5 for various M and X are shown in Figure 2.3. The dotted curves in
Figure 2.3(a) and Figure 2.3(b) are respectively associated with the above parameters

M = 100 and X = 0.425.

2.2.2 Saturated Stress and Nonsaturated Stress

Let T and a give the first Piola—Kirchhoff and Cauchy stress tensors, respectively.
They are connected by

a' = grrT. (2.16)

Mechanical equilibrium is governed by the usual equations
DivT = 0 on fix, 4: diva = 0 on Sly, (2.17)

where 9X and fly denote the gel domain in the reference and deformed configuration,
respectively. We continue to let Vpoly and Vliq be the volume of polymer and fluid

respectively. In particular

Vpoly = / dVX. (2.18)
Qx

After mixing the local gel volume is J so that the overall gel volume is

24

 

 

 

 

    

 

 

 

 

. ‘. x =‘0.000
221'. . x = 0.225“
'. —X = 0.425
,, ----- X = 0.625
2‘ ---X = 0.825
1.8-

v V
1.6 ~
1.4-. ................

1.2 ----------------------------
1

0 0.2 0450.6 0.8 1

(a) M = 100 and X varying from 0.000 to 0.825

 

 

]) ----- M/[J=2O

2.4":v "'M/#=5O J
' —M/n = 100

 
 
 

 

 

  

 

2.2" ‘ M/# = 200 l
. M/n = 1000
2. .
V l
1.8' 3

1.6'
1.4'

 

 

 

 

 

(b) X = 0.425 and M varying from 20 to 1000

Figure 2.3: The free swelling stretch ratio C as given by (2.14) and (2.15) is a function
of M, X and C. The dependence of C on C is shown for various M and various X

25

'" “’I

V E J dVX. (2.19)
9x

The maximum overall gel volume is V

poly + Vliq corresponding to uptake of all fluid

into the gel. So long as this does not occur the gel is able to saturate. The system is

not saturated only if V = Vpoly + Vh-q. Together this gives the global constraint

(J —- 1) de 3 V1,, (2.20)
Qx

If the system is saturated then the first Piola—Kirchhoff stress tensor is given by
T = 6W/8F and the saturated Cauchy stress tensor is
2 84>

=_ _ T
0' JFBCF +h(J)I, (saturated), (2.21)

where the first term on the right hand side follows exactly as from conventional
hyperelasticity and the second term on the right hand side follows upon recalling
that oJ/or = JF-T.

If, however, the material is not saturated, then a constant reaction stress —pI is
generated so as to enforce the constraint (2.20). This gives the nonsaturated Cauchy

stress tensor

0 = %F3%FT + h(J)I — pI, (not saturated). (2.22)

Thus, as in incompressible hyperelasticity, there is a pressure contribution to the
Cauchy stress tensor that is not determined by the deformation.

Here, however, it is important to emphasize that the constraint (2.20) is a global
one. This is in contrast to the situation in conventional incompressible hyperelas-
ticity in which the constraint detF = 1 is a local one. The pointwise constraint in

incompressible hyperelastic gives rise to a reactive pressure that can vary spatially. In

26

contrast, the pressure p in (2.22) takes the same value at every point in the nonsatu-
rated gel. In particular, since p in (2.22) is constant, a condition of spatially varying
J does not generally permit the last two terms in (2.22) to be consolidated into a
single term with redefined p.

Evaluating (NI/BC in terms of the dependence on 11,12, J gives

_ 26¢ 2 2 64> 6(1) 66
0—-751—2B +J(a—Il +116_I2)B+ (a—j‘f' h(J))I, (saturated), (2.23)

while the nonsaturated Cauchy stress tensor becomes

__ 2 04> 2 04> 84> 6CD
a— J012B +J(6—1—1+116—12)B+(E—J+h(J)-pl) , (notsaturated).-
(2.24)

For the uniform expansion F = J 1/3I, the saturated Cauchy stress (2.23) is a multiple
of the identity tensor, say a = —;13(J)I. Here 13(J) gives the mechanical pressure
associated with a given value of J and hence a given local fluid concentration within
the gel. The condition 13(st) = 0 recovers the free swelling condition (2.12).

A simple but fundamental application of (2.23) is to determine the infinitesimal
shear modulus G and bulk modulus K for a saturated gel. First consider a simple
shearing deformation after free swelling: y1 = CX1 + nCXg, y2 = CX2, y3 = CX3,
where K. is the amount of shear. For this deformation, the shear stress 012 for a

saturated gel is obtained from (2.23) in the form

28<I> 0(1)

012=G(n)n, C(K. )=— (59—1—1 + 2C—— (91.2 (2.25)

The infinitesimal shear modulus for a saturated gel is therefore given by

1 66 66
0(0) = 2 (Ea—I1 + (8—12) (2.26)

 

11=3C2,12=3C4.J=C3

27

Another important elastic parameter for saturated gel is the infinitesimal bulk mod-
ulus K. Consider a uniform compression: y = fiCX, where 8 is the compression
stretch after free swelling. The associated uniform compression stress a derived from
(2.23) is

2 8<I> 8<I> 8(1)

0 = 3287+ +4s4— + '67 + h(J). (2.27)

By definition bulk modulus is given by

K=—V—a—0—=—-83C3V 80 880

av (la—(3353(4)) = 755- (2-28)

Incorporating (2.23) into (2.28), and assuming that 11, 12, and J are mutually inde-
pendent in the expression of (I) (which is not unusual), the infinitesimal bulk modulus

for the saturated gel is then given by

264 4CB<I> C3[62——¢+ +d_h_(J)]

 

43— ‘ 34611 318.12 (M
2 (2.29)
56 (D
"K's—12771“ 612
2:11—342, 12: —,344 13:43

28

‘-

Chapter 3

Homogeneous Deformation

In this chapter we consider the homogeneous deformation response using the
constitutive forms (2.13) and (2.5) in (2.11). Since J is a constant for homogeneous
deformation, condition (2.20) reduces to the same condition as in free swelling, namely
(2.9). This in turn is equivalent to the condition J S J... where J... is defined in (2.10).

The Cauchy stress tensor then follows from (2.23) and (2.24) as
_ fl] 2
a _ j (1 — 4 +411)B — {B ] + h(J)I, (1 g J < J...), (3.1)

0’:

kl:

[(1 — E + 41013 — 4132] + (W...) — nu, (J = 1+). (3.2)

Since the term (h(J...) — p) in (3.2) involves constant h(J...) and an arbitrary constant
p, it follows that h(J...) can be absorbed into p. However in later chapters, when
nonhomogeneous deformation is considered, it will be necessary to retain h(J) in the
nonsaturated Cauchy stress tensor.

We now consider in turn: equibiaxial loading, uniaxial loading, and equitriaxial

loading. Each is referred to a coordinate system with mutually orthogonal unit vectors

e1, 62, e3.

29

3.1 Equibiaxial Stress

This is the specialization 0’ = o(e1 8) el + e2 <8 e2) and we restrict attention to
deformations that preserve the loading symmetry so that F = A(e1 (8) el + e2 ® 92) +

(J / A2)e3 ® e3. For a saturated material it follows from (3.1) that

p [2; +6 (g (A2 — 1) + 1%)] + h(J) = a, (saturated), (3.3)
and
[4 [Al4 + C (:3 (2A2 — 0)] + h(J) = 0, (saturated). (3.4)

The second equation provides the relation between J and A whereupon the first
equation provides a in terms of either A or J. The graph of stress vs. stretch passes
through the point (A,o) = (C, 0) which corresponds to the state of free swelling.
Recalling the requirement that J 2 1 consider the limit J —> 1 in (3.3) and (3.4). For
h in (3.4) given by (2.8) this limit gives

[—p(1i_51n(1-J_1)]—1/4, 1f0g4<1,
A ~ (3.5)
(-3411. (1 — J-1)]'1/2, if€= 1,

so that, in all cases, J —-> 1 as A —-> 0. According to (3.3) this requires a —+ —00 as

follows:

-u(1 -01—4, ifo s 4 < 1,

0 ~ h(J) ~ Mln (1 — J—l) ~ (3.6)

—2,n).—2, ifC = 1.

Thus the limit a —> —oo forces all of the liquid out of the gel, and in this limit A —+ 0.

30

Returning to the overall stress behavior, we find for M = 100p and X = 0.425
that J is monotone in A for all 0 S C S 1. The graph of 0 vs. A is also found to be
monotone for these parameters (Figure 3.1). Thus increasing stretch corresponds to
both increasing stress and increasing fluid content. This situation holds so long as
J < J...

A transition from a state of saturation to a state of nonsaturation occurs when

a:

e qbi’ and a umque value of stress,

J = J... This occurs at a unique value of A, say A
say Uqui' For A Z qubi the value of J remains fixed at J... and the relation between
stress and stretch is now determined on the basis of (3.2). This gives F = A(e1 <8)

él + e2 (8) e2) + (J.../A2)e3 8) e3 and

2 4
o = p [(1 —— C) (3— - 3%) + C (%_ —- 3%)] , (not saturated). (3.7)

It is to be remarked that the nonsaturated response does not depend upon h(J).
The solid curve in Figure 3.2 shows the saturated backbone curve 0' vs. A for
C = 0. The transition stretch qubz. can take on any value as determined by J..., and
three values of )‘qui are taken as examples. At each such Aquz. there is a transition
from saturated to nonsaturated response, and the associated nonsaturated response
curve that branches off of the backbone curve does so in a continuous but nonsmooth
fashion. The nonsaturated response involves an abrupt stiffening as evidenced by
the greater slope of the nonsaturated curve at the branch point. To be specific,
suppose A Z Aquz. so that the gel is not saturated because insufficient liquid is
available. Consider the difference in the value of stress for the nonsaturated response

and the value of stress for a hypothetical saturated response had sufficient liquid been

31

 

 

 

 

 

"28.5 1 1:5 2A2.5 a 3.5 4

(a) 47/14 versus A for equibiaxial loading

 

 

    
   

 

 

  

 

12 ' ‘ 1
_€ = 0.0
10‘ “f = 0.2
...... 5 = 0.4
8- IIIII 6 = 0-6
*6 = 0-8 ’H
5 61-62120 .....
4 — ---.:::::::‘.‘;‘.“.‘ :12 ........
2_ ..

 

 

 

8.5 1 1.5 g 25 a 3.5
(b) J versus A for equibiaxial loading

Figure 3.1: Stress-stretch behavior and volume-stretch behavior for equibiaxial load-
ing of a saturated gel with M = 100 and X = 0.425 showing the effect of different

32

 

 

 

1 _
3 °“
b -1.
-—satu'rated
-2 - - - - not saturated ‘

 

 

 

 

 

 

0.5 i 115 2 2:5A8 3:5 4 4:5 5

Figure 3.2: Stress-stretch behavior for equibiaxial loading with M = 100, X = 0.425,
and C=0. Three nonsaturated response curves are depicted. Each nonsaturated
response branches off of the saturated response “backbone curve”.

available. This difference is

an... - as... = [#(1- 042(71; —- })] + [Egg-90 — 1.)]

+[nga4(71; _. 9] + [£41 — .10],

where J in the above expression is the saturation value. Since J > J... if A > A211”.

(3.8)

it follows that each of the four separate bracketted terms in the above expression is
positive after loss of saturation. One therefore formally verifies that anon > Usat if
A > A2qbi' In addition, by differentiating the expression in (3.8) with respect to A and
then evaluating the result at A = qubz- and J = J..., one obtains the slope difference

at the location where a nonsaturated response curve branches off of the saturation

33

response curve. This calculation gives

d

2 4
d—A- (Unon’asat) [AZAak = #[(1—€)(i2—+;14)+€(3_3+11§)[[AzquM%[A=A:q

eqbi bi

(3.9)
Noting that the bracketted term in the above expression is positive, it follows that
the sign of the above expression is determined by the Sign of the derivative dJ/dA.
Since, as remarked above, J increases with A it follows that the slope difference is
strictly positive, thus quantifying the abrupt stiffening.

Returning to Figure 3.2 it is observed that two of the three nonsaturated response
curves depicted in this figure involve values of J... that give 47qu > 0. For these two
J... there is sufficient liquid for free swelling and the resulting loss of saturation occurs
once 0 increases to (7qu > 0. The third nonsaturated response curve depicted in
Figure 3.2 involves a value J... that gives (7qu < 0. In this case there is insufficient
liquid for free swelling to proceed to its saturated value. Thus the value a = 0 occurs
on a nonsaturated response curve. A biaxial compressive stress will, in this case, give
a response that proceeds down the nonsaturated curve until it joins the saturated
backbone curve. The resulting transition from nonsaturated to saturated response

now involves an abrupt softening in the Cauchy stress response.

3.2 Uniaxial Stress

This is the specialization a' = oel 03> el and we again restrict attention to defor-
mations that preserve the loading symmetry so that F = Ael (8') el + ./ J / A(e2 (8) e2 +

e3 8) e3). For a saturated gel it follows from (3.1) that

2

,u [(1 — C)? + 2CA] + h(J) = a, (saturated), (3.10)

34

 

and

u[(1-—C)31\-+ C(A + 31-15)] + h(J) = 0, (saturated). (3.11)

The latter provides the kinematic relation between J and A (Figure 3.3(b)) while the
former then gives the relation between uniaxial stress a and stretch ratio A (Figure
3.3(a)). The graph of stress vs. stretch again passes through the free swelling point
(A,o) = (C ,0). The predicted behavior associated with the possibility of squeezing
out all of the liquid in uniaxial stress again requires the consideration of J ——2 1. For
h given by (2.8) it is again found that A —» 0 as J —» 1. In particular, one finds that
A ~ 0([1n((J — 1)-1))-1/2)110 < 4 g 1, while A ~ 0([ln ((J — 1)-1))-1) if 4 = 0.
In all cases or —+ —-00 as A —+ 0. As regards the overall stress behavior, we find that o
is monotonically increasing with A when M = 100p and X = 0.425 for all 0 S C S 1.

For C = 0 it is found that the volume increases monotonically with axial stretch
(viz. Figure 3.3(b)). Thus, as in the case of equibiaxial stress, there will be a loss of
saturation when J = J... The stress vs. stretch response must then be determined
using (3.2). This gives

A2 1

J*
o = 14(7): — -/\-) (1 — C + C7), (not saturated). (3.12)

Even though we are currently considering the case C = 0 we have, for the sake of later
discussion, given the more general expression in the above equation. Returning again
to the case C = 0, the associated nonsaturated response curve can again be regarded
as branching off of the backbone saturated response curve. As shown in Figure 3.4,
each nonsaturated response curve is well separated from the backbone curve. Indeed

as A ——> 00 the C = 0 specialization of (3.10) and (3.11) give

1

J ~ (MG — X))1/2A1/2 + -3—- + 0071/2), (C = 0, saturated),

(3.13)

35

 

 

 

 

 

 

 

 

10

 

(a) 47/11 versus A

 

 

 

 

 

1G .
_._§=0.0
—-—4=0.2
----- 4:0.4

8” ------- C=0.6
---C=0.8
—4=1.0

’3 6

 

 

 

 

(b) J versus A

Figure 3.3: Stress-stretch behavior and volume-stretch behavior for uniaxial loading
of a saturated gel with M = 100 and X = 0.425 showing the effect of different C.

36

and

 

 

#1 3/2 _ 2M 1/2 _
o N _ 1 2A _ _ 2A + 0(A ), (C — 0, saturated).
(Mg ‘20) / 3M(1 2x)
(3.14)

Specifically, the Cauchy stress for the saturated gel is 0(A3/2) as A —> 00. In contrast,
(3.12) indicates that the Cauchy stress for the nonsaturated gel is 0(A2), which
accounts for the increasing separation between the response curves in Figure 3.4.

It is useful to note that the separation between the saturated and nonsaturated
stress response for the C = 0 gel is not present in the corresponding first Piola-
Kirchhoff stress. Recalling (2.16) it follows that this stress is T = Tel <8) e1 with
T = oJ / A. It then follows from (3.10) and (3.11) that the large stretch Piola—Kirchhoff

stress response of the saturated C = 0 gel is

)1/2A_3/2 + 00—2), (g = 0, saturated). (3-15)

TNflA-#(M(% -x)
Turning to the first Piola—Kirchhoff stress response of a C = 0 gel that is not saturated,

it is immediate from (3.12) that

J4

A2 , (C = 0, not saturated). (3.16)

T=,uA——/i

Thus the saturated and nonsaturated Piola-Kirchhoff stress response have the same
leading order behavior as A —> 00. Moreover, the separation between the saturated
and nonsaturated Piola-Kirchhoff stress response curves is 0(A“3/2) as the stretch
increases, and hence vanishingly small. Since first Piola—Kirchhoff stress directly scales
with applied load in a standard testing device, it follows that such a testing device may
have difliculty distinguishing between saturated response and nonsaturated response.

The situation is even more complicated for 0 < C S 1 since we find for this case

that J exhibits a local maximum, say Jmax, at a finite value of uniaxial stretch.

37

 

 

—saturated .’
4 - - - -not saturated .' ,

 

 

 

 

 

 

 

'20 0:51 1522; 8 354415 5

Figure 3.4: Stress-stretch behavior for uniaxial loading with M = 100, X = 0.425,
and C =0. Three nonsaturated response curves are depicted.

38

If Jmag; < J... then there is always sufficient liquid for the gel to remain saturated.
However if Jmaz > J... then the saturation value of J will be greater than J... on a finite

:1:
interval of stretch, say A um

—A < A < A* -B where the endpoint values A*

uni- uni— —A ’

AZm-_ B correspond to J = J... in uniaxial loading. On this interval the nonsaturated
response is given by (3.12). The nonsaturated response curve therefore branches off
of the backbone curve at )‘fmi— A and rejoins the backbone curve at )‘ftni—B' This
is shown in Figure 3.5, where it is to be remarked that the separation between the
C = 1 saturated and nonsaturated stress response is difficult to distinguish (unlike
the C = 0 case shown in Figure 3.4).

By a development parallel to that giving (3.8) one may show that any nonsat-
urated uniaxial response curve is above the backbone saturated uniaxial response
curve. The change in slope where a nonsaturated response curve connects to a satu-

rated response curve can also be found by the type of development that led to (3.9)

for the equibiaxial case. The corresponding result for uniaxial stress is

$920" ’030t)[1=1*:“[(1 —§)J_..+,\2il/\=A*n.dAl,\=/\;m' (3‘17)

uni

The Sign of this expression is again determined by the derivative term dJ/dA. This

derivative is positive at A*

um._ —A and negative at A*

um._ 8‘ Thus, under increasing
stretch, an abrupt stiffening occurs in the Cauchy stress for both the loss of saturation

transition at A*m_ —A and for the regain of saturation transition at Aum._ B'

3.3 Experimental Validation

In this section, we refer to R. Monroe’s experimental data obtained by uniaxial
stretch taets on hydrogels [74] in order to validate hyperelastic constitutive models

of the type presented here. Monroe constructed a uniaxial load frame that performs

39

 

 

 

 

 

 

 

 

3
2.9-
a 2.8-
2.7-
—saturated
- - -not saturated
2.6

 

1 1:5 2 2.5
A

(a) J versus A

 

   
   

 

—saturated
1'5' ---not saturated

 

 

 

 

 

 

 

 

 

 

0.5“ (anon _ asat)/Mq
E 0 x10'3 \
b 8 (”\l
-05. 6 "I \‘x
_. .. 4 '0' 'f\ “\‘
1 2 : '1' “\0‘ .‘0‘
_1.5 0 ”Amt”! ................ \
1.5 2
-2 L .
1 1.5 2 2 5

(b) o/u versus A

Figure 3.5: The dependence of volume and stress on uniaxial stretch for M = 100,
X = 0.425, and C = 1. Three nonsaturated examples are shown corresponding to loss

a * _
of saturation at )‘uni—A — 15,16 and 1.7.

40

individual or coupled experiments of stretching and swelling. The results shows the
monotone stress stretch behavior of hydrogel in the uniaxial loading tests (Figure
3.6(a)). Monroe found that the hydrogel volume increases with the applied load until
it achieves maximum swelling (Figure 3.6(b)). After that, further increase loading was
found to make the hydrogel shrink. This non-monotone volume change behavior is
consistent with the prediction of our hyperelastic constitutive model (c. f. Figure 3.3)
with a selected combination of material constants (0 < 6 S 1, M = 100, x = 0.425).

One of the main goals in [74] is to determine the material parameters from
both saturated and unsaturated testing, and then fit the stress-stretch constitutive
equations that are derived from the same free energy function as in (2.11). However,
to fit best with the experimental data, instead of using the form of (2.13), Monroe [74]

chose an alternative free energy model (Ogden Model) for the elastic free energy (I):
N p- a a a

<I>AAA= JAiAiAi- .1

(1.2.3)i__21ai(1+2+3 3) (38)

where Iii and a,- are material parameters obtained from fitting data.

In this work we use the Mooney-Rivlin form for the convenience of numerical
calculation, but it is the general approach to introduce different material models as
well as material parameters to make the theoretical prediction in consistent with the
realistic mechanical behavior of hydrogels that are composed of different polymers

and solvents.

3.4 Equitriaxial Stress

This is the specialization a' = 01 with F = J 1/ 31. The relation between a and
J then follows from (3.1) and (2.8) as

a = ,1[(1 — g)J‘1/3 + 2§J1/3] + M[J_1 + XJ_2 +1n(1— J—1)]. (3.19)

41

 

18000 . . .
16000
14000
12000

l—H
1

I
J.

I
1

10000 - 1
8000 ~ 1
6000 - I 1

0' [Pa]

4000 - ,
2000 ~ -
0 - . -

L

2 4 6 8 10 12 14 16
As

(a) Uniaxial stresses 0 versus saturated stretches As

 

 

 

 

100 ' ' F I I I I 1

90

80 l- l 1’
70 . 1 f 0 _

 

 

 

N.
60- ‘” J
50- -
40- I .
30

 

0 2000 4000 6000 8000 10000120001400016000
0’ [Pa]

(b) Volume ratios J versus uniaxial stresses 0

Figure 3.6: Stress-stretch behavior and volume-stress behavior for uniaxial testing of
a saturated gel as given by R. Monroe in [74].

42

Thus (3.19) with a = 0 and J = 1 /1/f3 retrieves (2.14). As one would anticipate,
a > 0 implies J > 1 /1/f3 and vice-versa. Consider again M = 100 and X = 0.425.
The graph of 0 vs. A = J 1/ 3 is then found to be monotonically increasing provided
5 > 0.038. However if 6 < 0.038 we find that the graph of 0 vs. A exhibits a stress
maximum (Figure 3.7). For 0 < § < 0.038 the stress maximum is followed by stress
minimum, after which the stress is once again monotonically increasing. For 6 = 0
the stress is monotonically decreasing to zero as A —+ 00.

Similar behavior occurs for values of M = M / ,a and X that are close to (M, X) =
(100, 0.425). Specifically, it is found that the graph of 0 vs. A = J1/3 is monotonically
increasing for all A only if g is greater than a critical value €cr2't- This critical value
varies with both A71 and X (see Figure 3.8). For certain values of M and X it is
found that {cw-it = 1 meaning that there always exists a local maximum in the stress
response for equitriaxial loading. This occurs for example if M = 100 and X 2 0.752
as can be seen from Figure 3.8(b).

As is well known in hyperelasticity in general, nonmonotone stress response be-
havior of the type encountered here when M = 100, X = 0.425, and 0 S 5 < 0.038 is
implicated in various loss of stability phenomena. Such issues are beyond the scope of
this dissertation. In contrast, for M = 100, X = 0.425, and f > 0.038 no such issues
arise. Then, J increases with a and loss of saturation occurs at the constraining value
J = J... when a = agqm. s u[(1 —§)J*_ 1/ 3+2§J3/ 3]+M[J:1+XJ;2+1n(1—J;1)].
For a Z azqtm. no additional expansion is possible, and the deformation gradient is
therefore stalled at F = J} / 31. The reactive pressure —p in (3.2) then renders the

equitriaxial stress 0 as formally indeterminant for J = J... (see dashed curves in Figure

3.9).

43

 

   

 

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

—g = 0.0 ‘
"-6 = 0.2 -
------ of = 0.4
-----f = 0.6 H
. 6 = 0.8 .
v g = 1.0
"8 1.5 2 2.5 3 3.5
A
(a) a/u versus A
0.8 . .
5 increasing from 0.00 to 0.05
with the increment
0-7'of0005
4
0.6r :
1 0.5-
\
b 0.4-
0.3’
0.2r
0.1 ‘ ‘
2 )1 4 6

(b) Magnified portion of (a)
F_igure 3.7: Stress-stretch behavior for equitriaxial loading of a saturated gel with

M = 100 and X = 0.425 showing the effect of different g. For E > 0.038 the graphs
are monotonic. For 0 < E < 0.038 the graphs are decreasing on a finite interval in 5.

44

 

0.1

0.08
30.00
“30.04-
0.02-

0

0

 

 

 

 

 

200 400 600 8001000
M / u

(a) {emit versus M with X = 0.425

 

0.8-

 

 

_

 

“0:2 04 0.6 0.8
X

(b) gcrz't versus X with M = 100

Figure 3.8: Critical value 5m.“ for the existence of a local stress maximum in the
equitriaxial stress response. For 6 > {crit the response is monotonic.

45

 

   
  
 

 

 

 

 

6* E : :
4* § 5 '
3. l
P 2
O .
—saturated
_2_ ---not saturated

 

 

 

1 115225 3 3.52145 5

Figure 3.9: Stress-stretch behavior for equitriaxial loading with M = 100, X = 0.425,
and (=05 showing three examples of nonsaturated response. Since the stretch A =
J1/3, no additional stretch is possible after loss of saturation.

46

3.5 Some Energy Considerations

We close this section on homogeneous deformation with some remarks upon
the energetic aspects of saturated response vs. nonsaturated response. A standard
expectation would be that a saturated response would minimize an appropriate energy
with respect to all nonsaturated responses, since the former is able to take on any
value of J. We consider this issue in the context of uniaxial response, since that case
offered the intriguing behavior such that, for certain 5, there could be a load induced
loss-of—saturation followed later by a load induced regain-of-saturation. Following the
previous development in Section 4.2, we again restrict attention to deformations that
preserve uniaxial symmetry and we also assume that there is sufficient liquid for a
saturated free swelling.

Consider a unit cube that is oriented with respect to (e1, e2, e3) prior to swelling
and prior to the application of load. After free swelling, let the two sides with normal
e1 be subject to a uniform normal traction with resultant normal force P and let
the four remaining sides be traction free. The cube is deformed into a rectangular
parallelepiped with stretch ratios A1 = A and A2 = A3 = m. One may distinguish
between either a hard loading device or a soft loading device. In the hard device, A
is regarded as given. For the soft loading device, P is regarded as given.

The energy to be minimized for the hard device is simply the stored energy
density as given by the integral of W(F) over the reference configuration. Since the

reference configuration is the unit cube, it follows from (2.11) that
Ehard = EhardO‘r J) E 009 + 2J/A, 2AJ + J2/A2, J) + H(J). (3.20)
For (I) given by (2.13) the expression for Ehard is

2
Ehardu, J) = §[(1—g)()12 + 3; —- 3) + €(2JA + $ — 3)] + H(J). (3.21)

47

For the hard device, one minimizes EhardO‘v J) with respect to J at fixed A. Formally,

the first step in the minimization involves seeking solutions to

3 _
37Ehard : 0! (3-22)

which gives an equation for J provided that J obeys J S J*. Notice that (3.22) with
(3.21) retrieves (3.11).

The soft device energy expression differs from the hard device expression by
subtracting the work that is done on the gel by the loading device. This represents
the change in potential energy of a conservative loading device so that the system
under consideration involves both gel and loading device. The work of the load P on

the unit cube is P(A — 1) so that
Esoft = Esoft(’\v J) P) E Ehard(’\a J) — P(A — 1)- (323)

For the soft device, one minimizes E30 ft(A, J, P) with respect to both J and A at
fixed P. Formally, the first step in the minimization involves seeking solutions to

(9 6

557301150, and aEsoft=0 (3.24)

Using E30 ft = Ehard — P(A — 1) it follows that (3.24) is equivalent to

a - 3 -
57Ehard = 01 and aEhard " P = 0' (3’25)

Since P = 0A2A3 = 0.] / A one verifies for (I) given by the Mooney-Rivlin form (2.13)
that (3.25) is the same as (3.10) and (3.11).
Thus both the hard and soft device analysis retrieves the saturation relation

(3.11) between A and J for uniaxial stress. This is a standard type of argument

48

in hyperelasticity, although in the present context of a saturated gel it permits the
determination of the fluid fraction (1 — 1 / J). In addition, the soft device analysis
gives the load that is necessary to support this deformation. These relations hold so
long as J S J... If, however, the solution to the above equations require that J > J...
then the minimization is found to be given by taking J = J... In this case the gel is
not saturated. Viewed in this fashion, the nonsaturated minimization occurs on the
boundary of the minimization domain. It follows that a nonsaturated solution can
never provide a lower minimum than that of the saturated solution, since the saturated
minimization is able to explore a larger space of J values. Such considerations arise
generally if a hyperlastic material is composed of a mixture of substances where one
or more of them (here the liquid) is in limited supply.

For specified A, let the solution to (3.22) be J = Jsat(A). Figure 3.10(a) displays
Ehard(A,Jsat(A)) vs. A for the standard constitutive forms (2.5) and (2.13) with
g = 1. This is the hard device energy associated with liquid saturation in uniaxial
stress. This panel also displays the nonsaturated hard device energies Ehard(’\1 J...)
that are associated with the transitions previously displayed in Figure 3.5. The
nonsaturated gels involve more energy than the saturated gel, although it is a small
effect in the figures. The saturated hard device energy EhardO‘a Jsat(A)) takes on its
smallest value at A = C which corresponds to free swelling.

In this context it is instructive to consider the graph of uniaxial load P as a func-
tion of A for both saturated and nonsaturated response, again using the constitutive
forms (2.5) and (2.13). The relation between P and A for saturated response follows
from (3.10) and (3.11) using P = aJ/A. The graph of this relation now supplies
a backbone curve. The corresponding relation between P and A for nonsaturated
response follows from (3.12) using P = 0J1. / A. Since the reference area of the unit
cube is unity, it follows that the load P is simply the first Piola—Kirchhoff stress T,

which we recall appeared previously in (3.15) and (3.16) for the case of the 5 = 0 gel.

49

 

 

 

 

 

 

 

 

°‘
- - -n0t saturated

 

  

 

 

 

 

 

 

 

 

-021 1.5 2 2.5
A
(a) Ehard versus A
3 .
—saturated
2 . - - -not saturated
1 . _
\ (Pnon _ Panza/Ill
o. _3 .
i x 10
\ ..1 n.‘
DH 5 I" “‘
—2 - :' ‘3‘
I ,' A
_3. O ‘31-“! "I
-4 -5 ‘\/’
1: 1.5 . 2
J
1 1 5 2 2 5
A

(b) P versus A

Figure 3.10: Comparison of a saturated and a nonsaturated gel under uniaxial load.
On the left is a comparison of the hard device energy expression. On the right is the
resultant force as a. function of stretch. As in Figure 3.5 the constitutive parameters
are given by M = 100, X = 0.425 and g = 1. The transitions are taken to occur at

50

For general 5 it follows on the basis of the above considerations that

Pm — 3a. =#(J—J*)[(1—o$+€(%l- —1)]. (326)

Here J = J (A) follows from (3.11), and so is given by the function J = Jsat(A) as
defined above.

Consider again the case that is associated with Figure 3.10 in which the nonsat-

*

urated response is confined to the interval A; m-_ A S A g A; ni— B where Auni— A

and Aim; B depend upon J... The graph of the P vs. A saturated backbone curve is
depicted in Figure 3.10(b) along with the nonsaturated P vs. A curves on the finite

intervals of nonsaturation. The nonsaturated curves depart and rejoin the backbone

*

uni—B' It is noted however that there is an intermediate

curve at )‘flni—A and A
value of A at which each nonsaturated curve crosses the backbone saturation curve.
The nonsaturated response is first above the saturated response and then below the
saturated response. It follows from (3.26) that the crossing value of A is a root to
(1 — £)/A2 + €((J... + J)/A3 — 1) = 0. This crossing value of A depends on J... and we
find that it occurs before the value of A that gives Jmax for saturated response.

Here it is useful to note that the ordinary derivative of Ehardfl, Jsat(A) with
respect to A gives, after use of the chain rule and application of both results from
(3.25), that

P = diAEhardO‘: J 3at(A)), (saturated). (3.27)

The corresponding result for the nonsaturated force-stretch response at any fixed

value J... is given by

P = deEhardmv J...), (not saturated). (3.28)

51

Combining (3.27) and (3.28) gives a generalization of (3.26) in the form

(I — _ _
Pnon - sat : a (Ehard(’\v J*) _ Ehard(A: Jsat(/\)))- (329)

We may use (3.29) to show that the crossing exhibited in Figure 3.10(b) is not
tied to the particular constitutive forms (2.5) and (2.13). Consider any H (J) and
<I>(I1,12,J) such that the saturated uniaxial response gives a local maximum in
the relation between J and A. In such a case, if J... as given in (2.10) is some-
what less than the maximum value of J, then there will be an interval of nonsat-
uration which will again be denoted by Azm; A S A S AZRi—B' In particular,

Jsat(A;m-_A) = Jsat(A;m-_B) = J... Integration of (3.29) now gives

 

Afini—B A*

(Pnon — sat) d’\ = (EhardO‘: J*) _ Ehard(’\1jsat(’\)) Aunt—B = 0 _ 0 = 0'
* uni—A
uni—A

(3.30)
It thus follows that there is zero net area between the saturated force curve and each
nonsaturated force curve on the finite interval of nonsaturation. Thus crossing of the
saturated and nonsaturated force-stretch response curves must occur within the finite

interval of nonsaturation.

3.6 Summary

In this chapter, we discuss homogeneous deformation of swollen elastomeric gels
based on the hyperelastic constitutive theory described in Chapter 2. To do this,
we distinguish between the liquid saturation and the situation of loss of liquid satu-
ration. Equibiaxial loading, uniaxial loading, and equitriaxial loading problems are

then considered in turn.

52

o For equibiaxial loading, it is found that for M = 100u and X = 0.425 that the
both 0‘ and J is monotone in A for all 0 S 5 S 1. Thus increasing stretch leads to
both increasing stress and fluid absorption as long as the gel is saturated. When
a transition from saturation to loss of saturation occurs, the gel’s overall volume
is then fixed (J = J...). Further increasing stress could continue monotonically
increase the equibiaxial stretch but resulting in a stiffer response than in the
saturated case. It is also found by analytical derivations that for the possibility
of squeezing all of the fluid out of the gel (J = 1), the equibiaxial stretch A has

to approach 0 and the corresponding stress much be negative infinite.

o For uniaxial loading, it is found that a is monotonically increasing with A when
M = 100p and X = 0.425 for all 0 S E S 1. The volume also increases
monotonically with axial stretch for 5 = 0. However for the cases with 0 <
g S 1 we find that the volume change ratio J exhibits a local maximum at
a finite value of uniaxial stretch. Corresponding experimental results are also
compared and qualitative consistence is obtained with the theoretical prediction
in this work. The non-monotone behavior could result in a phenomenon of the
nonsaturation bifurcation where the nonsaturated response curve first branches
off of the saturated curve and later rejoins it. Analytical derivatives again shows

that for all casesa——> —00 and A—->0as J-* 1.

o For equitriaxial loading (J = A3), the stress is found to be monotonically in-
creasing if E > Ecrit: and exhibits a stress maximum if 6 < {crit- This critical
value varies with both material parameters M and X. The bifurcation of stress
stretch curves is also discussed but it is trivial since no further equitriaxial

stretch is possible as the gel’s overall volume is fixed (J = A3 = J...).

o The attempt of explaining the load induced loss-of-saturation followed later

by a load induced regain-of-saturation in the case of uniaxial loading from the

53

energetic viewpoint is made. The energy-stretch relation is consistent with the

stress-stretch relation previously obtained.

54

Ii

Chapter 4

Eversion of a Swollen Cylindrical
Tube

The eversion problem of a hyperelastic tube was first presented by Rivlin [75]
for what is now known as the Mooney-Rivlin constitutive law. It was found that,
in general, zero-traction boundary conditions could not be achieved pointwise at the
two ends, for the assumed cylindrical deformation. This motivated Rivlin to employ
zero end resultant load conditions instead. Varga [76] gave experimental examples
for the eversion of different cylinders. Results showed that the the everted tube is
approximately cylindrical in the region away from the two ends.

Haughton and Orr studied the eversion of both incompressible [77] and compress-
ible [78—80] hyperelastic cylinders. The existence and uniqueness of the solutions to
such problems were discussed for the materials with a class of free-energy functions
proposed by Ogden [52]. In [78], an exact solution that satisfied zero—traction point-
wise boundary conditions both at the lateral surfaces and the two ends was obtained
for a particular compressible hyperelastic material. The specified strain-energy func-
tion led to a 1-D second order ordinary differential equation that could be solved
analytically. The analytical solution also shows an interesting result that the ratio
of the inner and outer deformed radii is equal to the ratio of the undeformed radii.
Several other material models including the class of ’Varga’ materials were considered

for this eversion problem in which case only numerical solutions could be obtained.

55

For the more general material models associated with the so-called semi-linear
strain—energy function there exists an analytical solution to the equilibrium equation
of cylindrical deformation. Specially, a new way of producing the analytical solutions
to this eversion problem was presented in [79]. This is involved replacing a number in
the equilibrium equation of the deformation r(R) by a constant parameter in order
to generalize the equation involving r(R). Then by working backwards, a class of
energy functions that would generate analytical solutions was obtained.

Furthermore, Chen and Haughton [80] studied the eversion problem for com-
pressible cylinders with the following boundary conditions: free traction at the two
ends and outer surface; either a zero traction (cavity remains open) or displacement
(closed cavity) boundary condition is specified at inner surface. They proved math-
ematically in [80] that exact solutions to the equilibrium equation with cylindrical
eversion deformation exist if two conditions are satisfied, namely: (1) the cavity re-
mains open during the eversion; (2) the strain-energy function of the hyperelastic
material can be written as the sum of three minor functions each of which depends

only on one of the three principal stretches.

4.1 The Inhomogeneous Deformation of an Everted
Tube under Axial Load

In this Chapter, we consider the inhomogeneous deformation of an everted tube
under axial load based on the hyperelastic continuum theory described in Chapter
2. For inhomogeneous deformation, the spatially varying stress field must satisfy the
equations of equilibrium (2.17). For the constitutive forms (2.5) and (2.13) the Cauchy
stress tensor is given by either (2.23) or (2.24) depending on whether the material is
saturated or not saturated. The only difference between these two formulae is due to
the presence of p in the nonsaturated Cauchy stress (2.24). Recalling that this p does

not vary with location, it follows that div( -pI) = —Vp = 0. Thus the equilibrium

56

equations diva = 0 are the same whether the gel is saturated or not saturated.
The distinction in mechanical response is due to the requirement (2.20) which causes
the gel to obey the constraint V = Vpoly + Vliq (viz. (2.18), (2.19)) when it is no
longer saturated. In the formal mathematical procedure, this constraint is met by
the presence of p in the boundary conditions. These issues will be demonstrated for
the inhomogeneous deformation of an everted tube subject to axial load. We shall

continue to use the Flory-Huggins constitutive form (2.5) for H and the Mooney-

Rivlin constitutive form (2.13) for (I).

4.1.1 Formulation of the Everted Tube Problem

Consider a hollow circular cylinder using polar coordinates in the reference con-

figuration 9X with
RingRo, OS9<27r, OSZSL, (R0>R,->0). (4.1)

Let it then be immersed in the fluid. In the absence of external tractions the gel
undergoes free swelling. Provided that sufficient liquid is available for saturation, the
deformation of the freely swollen gel cylinder is described in polar coordinates (7‘, fl, 2)

as

f:(& 0:9, zzgz am

Here C is again the free swelling stretch ratio as defined in (2.15).
Next, the swollen cylinder is everted (turned inside-out). The resulting defor-
mation is described with respect to polar coordinates (r, 6, z) obeying r,- S r S r0,

0 S 0 < 27r by
rsz) 0=a z=-Afi, dz>m. as)
The value A z is the mechanical stretch ratio in the axial direction, and will be called

57

mechanical stretch for short.
Combining (4.2) and (4.3) gives the transformation from the unswollen state

(R, 9, Z) to the swollen everted state (r, 6, z):
r = f(CR) = r02). 0 = a = e. z = —A.z. (4.4)

with

riSrSro, 0S0<27r, —lSzS0, (4.5)

where A; = AZC > 0 and l = AzL. The presence of the minus sign in the expression
for z in (4.3) gives rise to the eversion, and we seek solutions such that f’ < 0 so that

r' < 0. In particular, it follows that

r(R,) = r0, r(Ro) = r,-. (4.6)

For now, A; is regarded as a free parameter. It will ultimately be associated with the
axial load on the ends 2 = —l and z = 0.

The deformation gradient tensor associated with the mapping (4.4) is F = r’er ®
9R + (r/R)e0 <8) e9 — Azez (2) e2 where er, e9, ez and eR, e9, e2 are respectively the
cylindrical polar unit basis vectors in the deformed and reference configurations, and

r’ = dr/dR. It then follows that

11 = .12 + (32,-)2 + A3, 12 = (iffy + (r’xz)2 + (%,\z)2, J = -%",\z.
(4.7)

For the constitutive relations (2.13) and (2.5) it follows from (2.23) that, for saturated

response, there are no shear stresses and that the normal stresses are given by

R

”0:, [<1 - 0 +41%)? +51% + 148.22 ”(ng +1“ (1+ Jell’
(4.8)

 

 

 

Orr:—

58

 

 

009 = U. R ————,[(1—€)+€A2+§r'2]+M[A::,+X(X§—7J)2+ln (1+A:T,)], (4.9)

_#R)\z

 

 

0” = [(1 _ 0+0” ”(5)0 +Ml,\: 13' +X(A::r’)2 +1“ (1+ 1%)]

(4.10)

For nonsaturated response a common constant p is to be subtracted from these normal
stress expressions.

With the above stresses, the equilibrium equations associated with the t9 and 2

directions are satisfied identically. The equilibrium equation associated with the r

direction is

[<1 -€> News] +0.7: .._.2... J}

-%[<1-€>+a%l[r’2- (02] =0-
(4.11)

 

 

53-12-21.

     

The lateral surfaces of the everted tube are taken to be free of external tractions so

that

arr(r,-) = O, (rm-(r0) = 0. (4.12)

Rearrangement of (4.11) yields a second order ODE for r(R) in the form

2g%§(r'R—r)+ [(1—§)+§,\2]["T—7{’;Z]+ (7_rr)F(R.~.~i.)
2r’ [(1—§)+§(—2+A 2)]+ Afz—F (R,,rr, Az)

 

II
T =—

(4.13)

59

where the function F (R, r, r' , A z) is defined by

 

F(R,r,r',AZ) = — —

1J2 [(1_§)+§(§2+A3)] M 1 2X
J . (70—

and J is given by (4.7)3. Thus (4.13) with (4.12) provides an apparently well posed
problem for r(R). The solution describes an everted cylinder of saturated material
so long as the condition (2.20) is met. The solution is dependent upon constitutive
parameters: u, M, X, 5, and upon the stretch parameter A3.

The solution to this problem will generally give a radially varying an. In partic-
ular, a condition of 022 being identically zero is not anticipated for any combination
of parameters. This is connected to the observed behavior in eversion type defor-
mations on rubbery materials mentioned at the beginning of this chapter. Namely,
some flaring out of the top and bottom of an everted cylinder are observed in the
absence of any tractions on the top and bottom surfaces. Thus the observed traction
free deformations are not in accord with the simple eversion deformation (4.3). This
suggests that tractions on the top and bottom surface caps would be necessary to
sustain the simple eversion deformation. The deformation constructed on the basis
of (4.13) with (4.12) is consistent with these expectations.

The normal stress 022 generates an overall axial force on the everted cylinder.

Normalizing this axial force by the original surface area of the caps gives

27r T0

—-—— razzdr. (4.15)
W(Rg - R12) Ti

Pcap E

It is to be noted that r,- and r0 defined in (4.6) are unknown before the solution is
obtained. For fixed material parameters, Pcap will depend upon A z. The force-stretch
relation Pcap vs. A z is analogous to the homogeneous deformation relation in uniaxial

stress between P and A discussed in Chapter 3. Recall that this is given in general by

60

 

(3.27), where the specialization to the Mooney-Rivlin constitutive form uses (3.21). In
particular, one may then compare the relation between P and A obtained from (3.27)
with (3.21), to the relation between Pcap and AZ obtained on the basis of (4.15)
after the solution of (4.12) and (4.13). The former gives the force-stretch relation
for the cylinder before it is everted so long as it remains saturated. The latter gives
the force-stretch relation for the cylinder after it is everted so long as it too remains

saturated.

4.2 Numerical Results for a Saturated System

We have integrated the differential equation (4.13) after suitable nondimension-
alization by means of a fourth order Runge-Kutta method using the same material
parameters: X = 0.425, M = M / u = 100 that were employed in the homogeneous
deformation study. A shooting procedure was used to meet the two point boundary
conditions. The unswollen configuration was taken to involve Ro/R, = 2. Calcula-
tions were performed for different values of 5 and hence different values of the free
swelling stretch 5. For each such 5, the problem was solved numerically for a sequence
of AZ, which it is recalled is the mechanical portion of the mechanical stretch ratio
that appears in (4.3).

Radial deformation of the everted cylinder is shown in Figure 4.1 (a) and (c) for
the respective cases 5 = O and 5 = 1. Here the horizontal axis represents dimen-
sionless radial coordinates R E R/ R, in the reference configuration while the vertical
axis gives the dimensionless radial coordinates 1“ E r/(5R,) in the deformed config-
uration. As A; increases, it is found that both the deformed inner and outer radii
decrease, which is consistent with the expectation of a transverse contraction under
axial extension. Figure 4.1 (b) and (d) shows the radial variation of J (normalized

by 53) which gives the gel’s local volume change and hence the amount of fluid ab-

61

sorbed into the cylinder. It is to be noticed in panels (b) and (d) that the horizontal
axis represents the deformed radial coordinate F and so the curves start and end at
different locations on the horizontal axis. For each value of mechanical stretch Az,
the volume change is relatively greater near the outer periphery, indicating that the
volume fraction of liquid increases with radius in the everted cylinder.

Graphs of the normal stresses are displayed in Figure 4.2 for both 5 = 0 and
5 = 1 at a variety of mechanical stretches A 3. Note that arr < 0 in the interior of the
everted cylinder. The hoop stress 066 is found to be compressive in the inner portion
of the cylinder and tensile in the outer portion (which is broadly consistent with the
notion that the liquid is “squeezed” toward the outer portion). The axial stress azz
also exhibits radial variation and correlates with stretch A z in the anticipated fashion.
The qualitative form for all of the curves is shown to be quite sensitive to the value
of 5.

The relation between A2 and Pam as defined in (4.15) is shown by the solid
curves in Figure 4.3(a) and Figure 4.3(c) for 5 = 0 and 5 = 1, respectively. These
are compared with the uniaxial loading (without eversion) as represented by the
dashed curves in the same panels. For the neo—Hookean type gel (5 = 0), the everted
cylinder undergoes greater elongation at a given load than does the uneverted cylinder.
However, the opposite occurs for 5 = 1. Panels (b) and (d) of Figure 4.3 show how the
total cylinder volume varies with mechanical stretch, again for the respective cases
5 = 0 and 5 = 1. Specifically, the total volume as given by (2.19) is normalized by
the free swelling volume (solid curve), and compared to the normalized volume of
the uneverted cylinder (dashed curve). The total volume of the uneverted cylinder is
determined by the value of J in the homogeneous deformation for uniaxial stress. This
value of J was previously displayed (vs. overall stretch) in Figure 3.3(b) using the
now standard constitutive parameters M = 100 and X = 0.425. Recall for uniaxial

stress that the relation between J and stretch was not monotonic for 5 > 0. Thus

62

 

 

 
 
 

 

3 f
—Az = 0.5 5‘
. — = 0.5
'"A2 = 1.0 - -z
25 ....... 5‘2 _ 2 0 1 2.5 : ---%Z = 10 l
_____ .. : ' '”"'Az = 2.0
2:.~ AZ -- 4.0 ‘ : E-I-‘Az = 40
~~~~~~ 2, "a . ‘
IL ~~~~ ”X .i'
....................... ,5 r'!
1'52. ........................... 1.5»
1 ' """"""""" 7 1-
0'5 1.5 2 0'50 1 2
R F
(a) 1“ versus R for 5 = 0 (b) J/53 versus F for 5 = 0
3 .
-S\z = 0.5 1'4 _/.\ 0 5
"'A2 = 1.0 . “z = '
. 1 1.3 --- _ ,
2.5 ....... A2 = 1.5 ....... gz — 1.0 ,
----- Az=2.0 1.2» 5.2-1.5 ,
""‘Az = 2.0 ,5
.. 1.1 ,-'
x :
H 1_ :0 xx
0 ,9-
0'51 105 2 0-70 1'
R F
(d) J/53 versus '7‘ for 5 = 1

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(c) 1‘ versus R for 5 = 1

Figure 4.1: Radial deformation (7" = r/(CR,)) and local volume change ratio as a
function of radial position. As in previous figures, M = 100 and X = 0.425. The

cylinder geometry is such that R, = 2R,.

63

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

versus F for 5 = 1

(b) Urr

0

(a) arr versus F for 5

 

 

 

_Az

 

 

 

.1.
11111
if
ii
.1;

 

 

 

 

 

 

A

—Az = 0.5

 

A

"'Az = 1.0

A

“'A2 = 2.0

A

""‘A2 = 4.0

 

 

 

 

 

(d) 099 versus F for 5 = 1

=0

Ffor5

(c) 099 versus

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3 .
—Az = 0 5 . :21? = 0'5
6 "'Az = 10 ...... dz = 1'0
3‘. ....... A2 = 20 II! ;.-.§Z = 1.5
4 2‘. ..... AZ = 4.0 5 i" AZ = 2.0 <
1 i i'! :5 ’1'
\g 2 .......... \g :
b """""""" b If ,,,
0 ~ __________ .1 0* --------- "’1’
"'0 1 _ 2 3 ”0 1 _ 2 3
r r

(e) 0;; versus F for 5 = 0 (f) on versus F for 5 = 1

Figure 4.2: Radial distribution of normal stresses for the everted cylinder with R0 =
2R,- using the same material parameters as in Figure 4.1.
the dashed curve in Figure 4.3(d) is not monotonic. It is to be observed from panel
(b) of Figure 4.3 that the everted cylinder for 5 = 0 preserves the monotonic relation
between volume and stretch that was found in the uneverted cylinder. Similarly,
panel ((1) of Figure 4.3 shows that the everted cylinder for 5 = 1 gives a relation
between volume and stretch that is not monotonic as was the case for the uneverted
cylinder with 5 = 1. It is also interesting to note for both 5 = 0 and 5 = 1 that the
total volume of the everted cylinder is less than that of an uneverted cylinder at the
same value of stretch. Thus for both 5 = 0 and 5 = 1 the everted cylinder holds less
fluid than the corresponding uneverted cylinder.
The form of the volume vs. stretch curves in panels (b) and (d) of Figure 4.3
have immediate consequences as regards the possibility of loss of saturation. The
results mirror the corresponding results for the uneverted cylinder, since the everted

and uneverted curves in each panel display the same qualitative form. In general, the

65

 

1.5 . . .
—Eversion
1.4' ---Um'aa:ial Loading

 

h.)

 

 

Pcap/l’l

 

 

 

 

-1.'

 

 

 

-—Eversi0n 0_9. /
- - - U m'axial Loading r/
é 0'° 1 1:5 2
A}:
(b) V/st versus A; for 5 = 0

 

 

 

 

    

 

 

 

 

 

 

1 15
{\z
(a) Paap versus A; for 5 = 0
4 . . . 1.05
—Eversi0'n
- - - U m'aarial Loading
2 - .- __________________ ‘ 1
o ~
3. I ; "
\u 3“0.95»
8 .' §
DH _2_ II,
0.9- ‘
_4’ I. —Eversion
‘ - - -Unia:m'al Loading
i 1:5 2 2.5 0"" 1 15:5 2
A2: {‘2
(d) V/st versus A2 for 5 = 1

 

 

0.5

(c) Pmp versus 3‘2 for 5 = 1
Figure 4.3: Resultant axial force as a function of mechanical stretch and total volume
as a function of mechanical stretch for both the everted cylinder and the uneverted

cylinder using the same material parameters as in Figure 4.1. For the everted cylinder,

the reference geometry is again such that R0 = 2R4.

66

saturated solution holds provided that
«Azurg — r?) 3 «L023 — R?) + vhq, (4.16)

where it is recalled that A2 = 5‘25. We now consider differences that occur for the
respective cases 5 = 0 and 5 = 1 which stem from the different qualitative forms of
the curves in panels (b) and (d) of Figure 4.3.

Consider first the material with 5 = 0. The monotone curve for the everted cylin-
der in Figure 4.3b then indicates that an everted cylinder of the 5 = 0 material will be
saturated if ;\z is less than a critical value. Conversely, the everted cylinder will not
be saturated if 3‘; exceeds the critical value. The critical value of mechanical stretch
is dependent upon the amount of available fluid, and can be directly determined with
the aid of Figure 4.3(b). It is to be noted from this figure that, for a given amount of
fluid, the transition value of mechanical stretch for the everted cylinder exceeds the
transition mechanical stretch value for the uneverted cylinder.

Consider now the material with 5 = 1. Then the graphs of volume vs. stretch
for both the everted cylinder and the uneverted cylinder each have a local maximum.
Consequences for the uneverted cylinder were discussed in Section 1 and similar con-
siderations apply to the everted cylinder. Thus the everted cylinder will be saturated
for all values of stretch if the amount of available fluid exceeds that associated with
the value of the local maximum. If less than this amount of fluid is available, then
there will be an interval of stretch for which the everted cylinder is not saturated. It
is to be noted from Figure 4.3(d) that the value of the local maximum for the unev-
erted cylinder gives V > st where st is the free swelling volume of the cylinder
before eversion. Conversely, the value of the local maximum for the everted cylinder
gives V < st. Consider therefore two identical cylinders, each of which is removed

from its liquid bath after free swelling so that the amount of fluid is fixed at the free

 

1See the paragraph between equations (3.16) and (3.17).

67

swelling value. Let one of the two cylinders now suffer an eversion. This everted
cylinder will remain saturated for any value of axial stretch. In comparison, the un-
everted cylinder will immediately lose saturation if it is stretched a small amount.
However, if the stretch of the uneverted cylinder becomes sufficiently large, then it
too will regain saturation.

The above interpretation for a 5 = 1 gel is a consequence of the inequality
V < st holding for all values of stretch in the everted cylinder. More generally,
for 5 > O we find that the graph of stretch vs. overall volume continues to exhibit
a local maximum. However, as shown in Figure 4.4(a), the value of V at the local
maximum is found to be greater than the free swelling volume provided that this
positive 5 is sufliciently small. For these materials it follows that a quantity of liquid
just sufficient for free swelling gives rise to an everted cylinder that loses saturation
on a finite interval of stretch. The relation between force and mechanical stretch for
a saturated everted cylinder using the same values of 5 as in panel (a) of Fig. 4.4
is shown in panel (b) of this same figure. These serve as backbone curves for any
nonsaturated response. In order to display force vs. stretch curves associated with
loss of saturation it is necessary to solve the boundary value problem for an everted

cylinder that is not saturated. We consider this issue in the next section.

4.3 The Everted Cylinder that is Not Saturated

We now consider the consequences of a limited fluid supply for the eversion defor-
mation. Specifically, if the total volume of a saturated, everted cylinder is calculated
to be greater than Vpoly + Vh-q then there is insufficient fluid for the cylinder to be
saturated. The resulting nonsaturated solution obeys

liq

V
Mr?) —- r?) = (R3 — R?) + '5’ (4.17)

68

 

 

 

 

 

 

 

 

 

 

 

   

 

 

 

 

 

 

 

 

3 —5 = 0.0
e“; 0 ---g = 0.2 ‘
Q” ------- g = 0.4
5 ----- g = 0.6
. ---5 = 0.8
, *5, = 1.0
'50 2 i 6

(b) Pcap versus 5.3

Figure 4.4: Total volume after eversion as a function of mechanical stretch, and
resultant force as a function of mechanical stretch, for various values of 5, again using

M = 100, x = 0.425 and R0 = 2a,.

69

where it is recalled A2 = ;\z5. As discussed previously, the stress equations of equi-
librium are unchanged from that for a saturated system. Thus the displacement field
r(R) is still subject to (4.13) where F is given by (4.14). The normal stresses are
modified by subtracting the common constant p from the saturated stress relations as
previously given in (4.8)-(4.10). Eliminating this constant between the two conditions

in (4.12) now gives

 

 

 

_Rz'7"z" _ 3 2 2 — -Rz' Ra 2 Hi _
Azm [(1 a) + £( 12,-) + gig] + M[ 1213-71 + fibrin!) + 1n (1 + —Azr,-r,-')] _
R070, 7'_0 2 2 - ‘R0 R0 2 R0
_ A270 [(1 — E) + €(R0) + éAz] + M[—Az7‘o’rol 'l' X(W) +111 (1 + —AZT'0(’I‘0’)]),
4.18

where Ti, = 7" (12,) and ro’ = r’(Ro). The second order ODE (4.13) is to be solved
subject to (4.17) and (4.18), after which 19 is obtained from (4.12).

For values Vh-q that cause the saturated solution condition (4.16) to be violated,
we have solved (4.13) using a double shooting method so as to meet both (4.17) and
(4.18). The solutions are dependent upon the mechanical stretch 3‘3 and the value
of Vliq/rrL appearing in (4.17). Each such Vim/«L determines the range of 5.; for
which there is loss of saturation. As discussed in the previous section, the range of
nonsaturated 5‘; can be unbounded (as in the 5 = 0 example of Figure 4.3) or can be
confined to a finite interval (as in the 5 = 1 example of Figure 4.3).

Consider the material parameters associated with the 5 = 0 saturated solution
fields that were previously displayed in Figure 4.1 and Figure 4.2. The relation
between total volume and mechanical stretch for this material was displayed in panel
(b) of Figure 4.3. Since this relation is monotonically increasing, it follows that loss
of saturation could occur for any value of stretch. Following previous convention, this
transition value of stretch will be denoted )1: and the nonsaturated solution therefore

applies for 5.2 2 3;. Figure 4.5 displays the nonsaturated solution fields associated

7O

with this same material for the particular value of Vh-q that gives A: = 2.0. Thus
for this qu the curves in Figure 4.1 and Figure 4.2 corresponding to ;\z S 5‘: = 2.0
would continue to apply. However there would now be insufficient fluid to saturate
the everted cylinder when 5.; > 2.0 and so the ;\z = 4.0 curves in Figure 4.1 and
Figure 4.2 would no longer apply. Instead, the solution fields for 5.2 = 4.0 are as
depicted in Figure 4.5. The solution fields for the nonsaturated everted cylinder
for other values of mechanical stretch 5‘; > x); are also shown in Figure 4.5. The
nonsaturated solution is identical to the saturated solution at the transition stretch
51; as is verified by comparison of the curves for 5‘; = 2.0 in Figures 4.1, 4.2 and 4.5.
As the mechanical stretch increases, comparison of Figure 4.2 and Figure 4.5 reveals
similar qualitative trends in the nature of the stress fields after loss of saturation.
In contrast, comparison of Figure 4.1 and Figure 4.5 shows qualitative differences in
the J field after loss of saturation. Specifically, the saturated solution gives values
of J that increase at each point in the everted cylinder. Hence all locations absorb
additional fluid as )1; increases so long as the gel is saturated. However, upon loss
of saturation, it is found that the value of J now decreases on the outer portion of
the cylinder (while continuing to increase in the inner portion of the cylinder). This
is shown in detail in the first panel of Figure 4.6. Thus, after loss of saturation, the
fixed amount of fluid within the cylinder exhibits a quasi-static migration from the
outer portion of the everted cylinder to the inner portion of the everted cylinder under
increasing mechanical stretch :\z.

In a similar fashion Figure 4.7 exhibits nonsaturated solution fields for the 5 = 1
material whose saturated solution fields were previously exhibited in Figure 4.1 and
Figure 4.2. The relation between overall volume and mechanical stretch for this
material is displayed in panel ((1) of Figure 4.3. Since this relation exhibits a local
maximum it follows that loss of saturation occurs on a finite interval of stretch. Again,

following previous convention, we denote the endpoints of this interval by X:_ A and

71

 

 

 

 

 

 

 

 

 

 

 

 

 

  
 

 

 

 

 

 

2
1.8*
1.6-
12' _;\z = 2.0 ‘
"':\z = 3.0
1* """" 32 = 4.0 ‘
""" ;\z = 5.0
08.5 1 1.5 2

$

(b) J /53 versus 7“

72

 

-0.01i ‘.
-0.03’

-0.05*

Orr/IL"

-0.07'

 

 

 

 

0.5 1 1.5 2

 

 

 

 

 

 

 

 

 

'1 1 _1f5 2
7‘
(d) 099 versus 7"
15 r
110' ‘~.
\ ......
N .....
N "
b
5’ ..‘~
\
0 1 _1.5 2
7'

(e) azz versus F

Figure 4.5: Deformation and stress fields associated with loss of saturation for an
everted cylinder with R0 2 2B,. The material parameters are M = 100, X = 0.425,
and 5 = 0. Loss of saturation is taken to occur at x); = 2.0. Thus the curves for
:\z = 2.0 are the same as the corresponding curves for 5 = 0 that were displayed in
Figures 4.1 and 4.2. For 3.3 > 2.0 the curves in this figure differ from corresponding
curves for a saturated gel.

73

 

 

 

 

 

 

 

 

 

 

0-7 transition >
§/ L Az—B
0.5 S21, 1.5 2 2.5
A2
0)) 5 = 1

Figure 4.6: Local volume change J versus 5‘; for an everted cylinder with R0 = 2R,-
and material parameters M = 100, x = 0.425 and either 5 = 0 or 5 = 1. On the left,
5 = 0 and the loss of saturation corresponds to that depicted in Figure 4.5. After loss
of saturation, increasing stretch redistributes the finite amount of fluid from the outer
region to the inner region. On the right, 5 = 1 and the loss of saturation corresponds
to that depicted in Figure 4.7. For the nonsaturated cylinder, increasing stretch now
redistributes the finite amount of fluid from the inner region to the outer region.

74

32—3 For the purposes of Figure 4.7, we take A}; A = 0.8 which in turn makes
A:_ B = 1.993. Once again, the fixed amount of fluid redistributes itself in a quasi-
static fashion as 3.2 varies in the interval A2; A S 5.2 _<_ 51:; B' In particular, as AZ
increases through this interval we find that fluid migrates from the inner portion of
the everted cylinder to the outer portion. This is shown in detail in the second panel
of Figure 4.6. For this material, after a return to saturation at X:_B, any further
increase in mechanical stretch gives a loss in overall volume, indicating that increasing
mechanical stretch now expels fluid from the everted cylinder.

Recall for the fluid saturated everted cylinder that the relation between axial
force and axial stretch was given previously in Figure 4.3 for both the 5 = 0 material
(panel a) and the 5 = 1 material (panel c). Loss of saturation causes the response
to depart away from this backbone response. The question therefore arises as to
whether this departure takes place in a manner that is similar to that which occurs in
the homogeneous deformation of uniaxial stress. Recall for the homogeneous defor-
mation of uniaxial stress that the departure of the nonsaturated axial force response
from the saturated backbone response was displayed in Figure 3.10(b) for the 5 = 1
material. The conspicuous feature of Figure 3.10(b) was that the nonsaturated axial
load response for the 5 = 1 material was first above, and then below, the saturated
axial load response. We find that similar qualitative behavior occurs for the 5 = 1
material as the everted cylinder is mechanically stretched. This is shown in Figure
4.8 where the difference between satuated and nonsaturated response is shown not
only for the case in which X:_ A = 0.8, but also for two other cases corresponding to

A:_A = 0.9 and A:——A = 1.0.

75

 

 

2.5 ' ..
—,\z = 0.8
"-33 = 1.2
2\ ....... :\z = 1.6 1
""" :\z = 1.993

Is1.5~

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.
0.51 15 2
R
(a) 1" versus R
1.4 ' . T
13* .i ""‘z = 1'2 H
; ....... A2 = 1.6
1.2» ----- ,‘\z = 1.993“
,,, 1.1-
x 5
'5 1_
0.9'
as
0.70 1‘ __ é a

(b) J/(3 versus 7‘

76

 

 

 

 

 

 

 

 

 

0.5 1 1:5 2 2.5

 

 

 

 

10
i
i
i
5 i
» l
1 f :
\ .I I
N I" I,
N j I
b 0 ..:"~.‘ 'I’
-5

(e) 07,; versus 1‘

Figure 4.7: Deformation and stress fields associated with loss of saturation for an
everted cylinder with R0 = 211,-. The material parameters are M = 100, x = 0.425,
and 5 = 1. The value of Vh-q is such that loss of saturation occurs on the interval
osgizgrma

77

 

 

— saturated

2 "'not saturated f
\ (Pa? — Pfé‘fiflu
- /- . .

 

 

 

 

 

 

 

 

 

1 112 1.4 1.61:8
1.5 2 2.5
A2

Figure 4.8: Resultant axial force as a function of mechanical stretch for the everted
cylinder with R0 = 2R.,-. Material parameters correspond to those in Figure 4.3(c),
namely: A7! = 100, x = 0.425, and 5 = 1. Three different values of Vh-q are considered,
corresponding to A: = 0.8,0.9, and 1.0. The nonsaturated response curve is first
above, and then below, the saturated response curve, and so mirrors the case of this
material under homogeneous deformation (Figure 3.10(b)).

 

 

 

78

 

4.4 Summary

The problem of everting a swollen tube subject to axial loading is studied in this
chapter. The inhomogeneous deformation response of the axial loaded everted tube is
compared to the homogeneous deformation response of the axially loaded tube when

it is not everted. Both systems can transition from saturation to nonsaturation.

c As AZ increases, it is found that both the deformed inner and outer radii de-
crease, which is consistent with the expectation of a transverse contraction under
axial extension. For each value of 31;, the volume change is relatively greater
near the outer region, indicating that the volume fraction of liquid increases

with radius in the everted cylinder.

0 The overall volume is calculated for the everted tube, and compared with those
for an uneverted tube. It is observed that the everted cylinder for 5 = 0 pre-
serves the monotonic relation between volume and stretch that was found in
the uneverted cylinder, and similarly the everted cylinder for 5 = 1 gives a
relation between volume and stretch that is not monotonic as was the case for
the uneverted cylinder with 5 = 1. More generally, we found that the graph of

overall volume vs. stretch exhibits a local maximum as long as 5 > 0.

o The transition to loss of saturation is also considered for the axial loaded everted
tube. For 5 = 0 it is found that after loss of saturation the fixed amount of
fluid within the cylinder exhibits a quasi-static migration from the outer portion
of the everted cylinder to the inner portion as the axial stretch increases. For
5 = 1, since the volume-stretch relation exhibits a local maximum it follows that
loss of saturation occurs on a finite interval of stretch. As ;\z increses through
this interval we find that fluid migrates from the inner portion of the everted

cylinder to the outer portion.

79

Chapter 5

Twisting of a Cylindrical Tube

In this chapter we discuss the absorption and redistribution of fluid in hypere-
lastic gels due to twisting; and in particular distinguish between saturated and un-
saturated gel systems. To this end we consider a boundary value problem of radial
displacement combined with azimuthal shear for an annular cylinder consisting of
a fluid infused hyperelastic media. The effect of either boundary displacements or
boundary tractions is considered so as to study how this alters the uniform fluid dis-
tribution when the hollow cylinder is in contact with a fluid bath at both inner and
outer radii. Certain aspects of this problem were previously considered by Rajagopal
and Wineman [48] for the case in which the overall volume is held fixed. Here we con-
sider certain generalizations in which the lateral surfaces may undergo both prescribed
radial displacement and prescribed relative twist. In particular these generalization
permit volume change so as to allow for the consideration of both saturated states
and unsaturated states.

In order to compare with the numerical results in [48] it is convenient that the

elastic energy (I) is chosen to have the same form as in [48]:

<I>(11,12) = 5‘22 [11 — 3 + a(12 — 3)], (5.1)

80

where he, > 0 and a 2 0 are material parameters. The case associated with a = 0
corresponds to 5 = 0 and a = 2110 in (2.13), whereas the case associated with a = 1

in (5.1) can be retrieved by setting 5 2 1/2 and [1 = 2M: in (2.13).
5.1 Twisting of a Saturated Tube

Consider a hollow circular cylinder of the gel described in the previous section.
Employ polar coordinates (R, 9,Z ) in the reference configuration such that R,- ;<_
R g Ro,0 S 6 < 27r,0 3 Z s L with R,- > 0. This cylinder is immersed in a
liquid bath and, following Rajagopal and Wineman [48], let the coordinates (1", 6, 2)

describe the mapping associated with this free swelling:
f = (R, 0 = 9, :2 = (Z, (5.2)

where C = J}:3 is the free swelling stretch ratio. For H given in (2.5) and (I) in
(5.1), the free swelling volume change J f3 is determined by (2.4). Wineman and

Rajagopal [48] using (5.1) for various 0 take
M = 2.379E7 dyne/cmz, pa = 2.375E6 dyne/cmz, x = 0.425 (5.3)

for a vulcanized rubber—toluene mixture. Figure 5.1 shows the the free swelling volume
change J f3 and the free swelling stretch 5 as a function of a for the parameter values
in (5.3).

In what follows it is assumed that there is always sufficient liquid available to
support free swelling. In view of the conditions (2.18) - (2.20), this is a requirement
that

V1,, > (st — 1)V,,0,y = (<3 — 1)7r(123 — 12,2)L,

and the volume of the freely swollen cylinder is st = J f szoly = (37r(R3 — R22)L.

81

 

 

 

 

 

2.

 

C

 

 

1

 

0 02 0.4006 0.8 1

Figure 5.1: Free swelling volume change J f s and stretch 5 for different a on the basis
of (2.3) and (5.1) using (5.3).

82

Now consider loading of the swollen cylinder. Following Wineman and Rajagopal
[48] attention is restricted to axially symmetric deformations that are described in

polar coordinates (7', 0, z) obeying r 2 O, O S 9 < 271' by
T = f(r), 0 = 5+ 5(f), z = X22. (54)

Combining (5.2) and (5.4) gives the deformation from the unswollen state with coor-

dinates (R, 6, Z) to the swollen loaded state with coordinates (r, 0, 2) as follows:
7' = f((R) = r(R), 6 = (i + 3((R) = e + g(R), z = AZZ, (5.5)

with AZ = izc.

The deformation gradient tensor for (5.5) is given by F = r’er 8) e R + (rg’)e9 <8)
9R + (r/R)eg (8) e9 + Azez <8) eZ where er,eg, e; and eR, e9, e2 are the cylindrical
polar unit basis vectors in the deformed and reference configurations, and 1" = dr/dR,

g’ = dg/dR. It then follows that

, I
we, 12 = (yaw/152+(gizizm'mz. J -—- T552-

(5.6)

2

,2 7"

I:
17" R

+(7‘9') +(
These expressions for 11, 12 and 13 match those given in (23) of [48] once allowance
is made for the slightly different notations used here.

For the remainder of this section it is assumed, unless stated to the contrary,

that the cylinder is saturated. It follows from (2.21) that the shear stresses on and

092 vanish. For constitutive functions given by (2.3) and (5.1) the other stresses are

Ritz RAz 2 1%
7'7" +X(rr’) +ln(1— rr’ )j’

(5.7)

given by

 

 

 

Rr 7' 2

 

83

 

 

 

 

 

 

 

(5.8)
R 2 2
ozZ = Ila—72‘ [1 + ar'2 + a (rg’) + a (%) J
2
+M [Inf +x (RA?) +ln(1— 1222)],
7‘?“ rr 1'1"
(5.9)
R9’ 2
07.0 = [la—[C (1 + (1)12).
(5.10)

With the above stresses, the equilibrium equation (2.17)2 in the z direction is satisfied

automatically, while the equilibrium equations in the r and 0 directions give

2
RA; + X RA; +1n 1_ RAZ
7‘7" rr’ 7‘7"

<1 mg) [7, — (9)2 — at] = 0,

(5.11)

 

 

 

(1 RT" 7‘ 2 2
87-{110—7 [1+a(§) +01%] +M

 

 

 

II I I
— _ - 5.12

where g” = d2g/dR2. It is to be remarked that one could obtain equations identical
to (5.11) and (5.12) using the framework in [48]. Namely, substitution from (25) -

84

(29) of [48] into (30) and (31) of [48] can be shown to result in (5.11) and (5.12) as
given above once allowance is made for the different notation.

We comment that, in contrast to [48], the methodology used here does not invoke
the mixture theory framework which, among other things, introduces separate mass
densities for the intermingled fluid and solid components as well as an interactive body
force when there is relative motion between these components. We do not require use
of this more general framework because the equilibrium problems under consideration
here always involve a stationary fluid component so that there is no relative motion
between the liquid and polymer components within the gel. If we further consider
time varying boundary conditions in this context, then the solutions that we obtain
describe a family of equilibrium solutions. In this case the fluid redistributes itself
in a quasi-static manner. This is the sense in which the term redistribution is used
in [48] and is also the sense in which the term is used here. As discussed recently
by Back and Srinivasa in [71], the slow diffusion of a fluid through a swelling (and
saturated) solid is also amenable to such a treatment. In fact a specific form of the
governing equation is given in (42) of [71] as

61$

div (8—FFT + 51) = 0, (5.13)

where 1,5 is the specific Helmholtz free energy and is akin to
E = W + .ulz'q(J _1)+ “poly (5-14)

in the present treatment. Here W is again given by (2.11) and “liq and “poly are
positive materials constants that correspond to the chemical potential of fluid and

polymer components, respectively. However the free energy 1; in [71] is defined per

85

unit current volume so that the connection to our E is that

113 = E/J (5.15)
whereupon
at} _ 16E 1 _T
55 — 755 7 7’”

by virtue of BJ/BF = J F_T. Hence it follows that

811; T ~_13E
F +1/2I—J6F

T _ .

whereupon (5.13) is identical with (2.17)2 since ”liq is constant.

5.2 Twist at Constant Overall Volume

We first consider the type of boundary value problem studied by Wineman and
Rajagopal in [48], which is stated as follows: the hollow tube after free swelling is
held fixed at its inner radius, i.e. r(Ri) = (R,- and g(R,') = 0, while the outer radius
is then rotated through a twist angle «p = g(Ro) while maintaining the same radius
r(Ro) 2 CR0. The governing equations are given by (5.11) and (5.12). Note that the
latter of these is, in its original stress formulation,

$1212 + 20_T9 = 0, (5.16)
dr 7"

which in turn can be immediately integrated to

T

are : 27TT2 ’

(5.17)

86

where T is the constant of integration. Since 27rR220r9| R1 = 277R30T9|Ro = T it
follows that T is the associated twisting moment, or torque, on the lateral surfaces.

A useful nondimensionalized torque T* is defined by the relation

,, T

= —————. 5.18

The boundary conditions are now taken to match those in equations (34) and (35)
of [48] meaning: (A) there is no additional radial displacement of the lateral surfaces,

i. e.

7"(RD = CR1: r(Ra) = CR0; (5-19)

(B) that the inner surface is rigidly clamped, i. e. g(R,) = 0; and (C) that the
outer surface is twisted by an an amount #1, i.e. g(Ro) = 11;. In fact, as regards the

boundary conditions on g(R), it is only the relative twist that matters, i.e.

g(Ro) - 90%) = 1P, (520)

since the addition of an arbitrary constant to the function g(R) merely represents
a rigid body rotation about the z-axis. An alternative to specifying 1/1 is to specify
the twisting moment T. In either case, the torque-twist relation (T vs. 7,11) is then a
characteristic feature of the macroscopic response.

Finally, in keeping with [48], we specify that the height of the cylinder is kept
fixed at the free swelling value f = 5 L. Thus the axial stretch AZ is the same as the
free swelling stretch 5 so that 5.2 = 1 in (5.4). This will in general require a resultant

axial force at the caps z = O and z = f given by

7‘0
P = 277/ ozz(r)7‘dr. (5.21)
T

i

87

 

It is to be remarked that an alternative problem is to find AZ, the value of the
mechanical stretch ratio, so as to give a zero resultant force at the two caps. We,
however, do not consider this possibility here.

With the inner and outer radius fixed at their free swelling values and the height
of the cylinder unchanged, the overall volume of the gel remains the same as that of the
free swelling state. This therefore results in no net intake or exit of fluid. Rajagopal
and Wineman observe this in their remarks regarding impermeable barriers at r, and
r0 before equation (18) of [48].

With the material parameters (5.3), in a cylinder with R0 = 2R,-, we used a
shooting method to obtain r(R) and g(R) obeying field equations (5.11) and (5.12)
and boundary conditions (5.19) and (5.20). The numerical results for the fluid redis-
tribution, interpreted in terms of the gel’s local volume change J — 1, were found to
exactly match those given by Wineman and Rajagopal in [48]. Among the main fea-
tures obtained by Rajagopal and Wineman [48], and confirmed here, are those shown
in Figure 5.2 where we have introduced the dimensionless radial position R = R/Ri.
Figure 5.2a and Figure 5.2c show the radial variation of J / J f3 for a = 0 and a = 1
respectively. This ratio decreases near the inner radius and increases near the outer
radius, which means that fluid redistributes itself from the inner region to the outer
region of the cylinder. The rotation of the cylinder is shown in Figure 5.2b and F ig—
ure 5.2d taking g(R)) = O. The angle of rotation increases rapidly near the inner
radius and gradually near the outer radius. The other figures in [48] can be similarly
reproduced including Figure 3 and Figure 7 of [48] which give the ratio of the current
liquid density in the gel to the free swelling liquid density. This ratio is denoted by

pf/pgat in [48], and the formal connection to our notation follows with the aid of

88

 

equation (24) in [48] as

pf/pfa. = <1— J‘H/u — 5.1). (5.22)
(density ratio from [48])

The normal stress distribution on (r) is not displayed in [48] and so, in particular,

the boundary tractions orr(r,-) and orr(r0) needed to sustain this deformation were
not given there. In Figure 5.3 we show the arr stress distribution. For a non-zero
twist angle 11) we find that or,- is monotonically increasing such that Urr('ri) < 0 and

(rm-(7‘0) > 0. Thus the trend in an- is similar to the trend in J. More generally,

graphs of the other normal stresses 096(R) and 022(R) are shown in Figure 5.4.

5.3 Twist Induced Fluid Desorption

It follows from Figure 5.3 that if r,- is held fixed and a twist angle a is imposed
at 7' = r0 then, for To to be kept at its free swelling value, it is necessary to supply
a tensile normal stress arr at r = 1'0. If such a normal traction is not present, then
one would anticipate a lesser value of re on the basis of a tensile traction having been
released. This motivates the consideration of a modified boundary value problem in
which the inner radius is fixed and twist is again applied at the outer radius, but the
outer radius is now taken to be free of normal traction. Thus boundary conditions

(5.19) are replaced by
r(Ri) = CR1, 01-7-(7'0) = O. (5.23)

Quantities associated with the solution of this problem are shown in Figure 5.3 where
(a)-(h) graph the same quantities that were previously shown in Figures 1-8 of [48]
respectively so as to show the effect of the outer surface being free of normal traction.
In particular (b) and (f) of Figure 5.3 confirm that r(Ro) < 5R0 indicating that the

overall volume has indeed decreased from its free swelling value. It is to be remarked

89

 

2.5

 

—T* = 0.000
---T* = 1.000 I,
2 ------- T’“ = 2.000 _1‘

----- T“ = 3.000
1,5 - T* = 4.000 .
- T“ = 5.000 5; .........

.u-u
............
..-

     
 

 

 

 

 

 

 

A
:UI'O,’

(a) J/st vs. R for a = 0
1.05

 

 

 

 

'3“: i. :3 0'.
;°-95 _.-',.-' —T* = 0.000
-' -' ---T* = 1.000

 

 

 

 

 

 

' .- ------- T‘“ = 2.000

0-9'5' ----- T* = 3.000

5 - T* = 4.000

0 85’ - :1“ = 5.000
' 1 1; 2

R
(c) J/st vs. R for 0 =1

D)

 

 

 

 

 

’.

 

 

 

 

1 1:5 2

(b)g=0-Ofora=0

1
-T" = 0.0

---T’ = 1.0
0.8 ------"T = 2.0 1
""" T‘ = 3.0 /
_ - T" =40 ,
0'6 - T“ = 5.0

 

 

 

 

 

 

 

 

 

,’
.0"
.0"
.¢
,4
"""
0 4 ..
. I
.I‘ .............
.. ’0 ...........
I - .....
I '- "’
n.’ - 0.
0 2. :I..I.’ ....
' .°.-,- .........
‘I’ ..........
I. . ”o
v"
G I
1 1.5 2

(d)g=0—Ofora=1

Figure 5.2: J normalized by the free swelling value J f3 = (3 as a function of R =
R/R, (panels (a) and (c)); and the twist g = 9 — 9 as a function of R (panels (b)
and (d)), for boundary conditions given by (5.19) and (5.20) with Ro/R, = 2. The
problem is equivalent to that presented in [48] and panels (a), (b), (0) correspond

respectively to Figures 1, 4, 5 in [48].

 

 

 

 

 

 

 

 

 

 

 

O
.-----'-'-"'-.-
.-
.o

I
.........
I
a"

"

    

 

—T* = 0.000 ‘
5' ---T* = 1.000
------- T* = 2.000
-2-_-' ----- T* = 3.000 '
..' ' T* = 4.000

- T* = 5.000

1 125 2
R
(b) a = 1

 

 

 

 

 

 

Figure 5.3: Normal stress or,» for the same problem as described in Figure 5.2.

91

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(b) 022 vs. R for a = 0

92

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 L5 2
R
(d) 0;; vs. R for (1 =1

Figure 5.4: Normal stresses 099 and on for the same problem as described in Figure
5.2.

93

that an overall volume decrease was also found to be the case in the torsion problems
discussed in [49] and [45].
1.2 1.05

 

 

 

 

.............
...........
........
.....
,,,,,,,,,,,,
""""""""
.......
.....
_,
.4
'

   

 

J/st

   

—T‘ = 0.01 \
---T‘ :10
------- T' = 2.0
-----T‘ = 3.0 .
. T' = 4.0
- T' = 5.0

 

 

 

 

 

 

 

 

 

 

0"1 1.21.4_1761:8 2 '1 1.21.4_1:6118 2
R R

(a) J/st vs. R for a = 0 (b) 'f/R vs. R for a = 0

 

 

 

 

 

 

 

 

 

 

 

0'41 1.2 1.4j6 1.8 2 1 1.2 1.4_1.6 1.8 2
R R

(c)pf/p£atvs.Rfora=0 (d)g=9—9vs.Rfora=0

Figure 5.6(a) and Figure 5.6(c) show the decrease in 17st as 11” increases for the
two cases a = 0 and a = 1. The dashed horizontal line in these figures corresponds
to the volume of the unswollen polymer. It is to be noted that this horizontal line is
consistent with it being the asymptotic limit of the upper curve in the event that all

of the fluid is expelled as 1/1 —» 00. Figure 5.6(b) and 5.6(d) give the corresponding

94

 

 

 

 

 

 

 

 

 

 

1.05
1 a
‘ 10-995'1'8"‘~. ...................................... :
Q: ........... __________________
3 0.99 =\_/'
. \ ,
L '-
0.985»
0.85-" . '
t? 0.98 \/
0.8 . , . . . . . .
1.2 1.4_1.6 1.8 2 0'9751 1.2 1.4_1.6 1.8 2
R
(f) f/R vs. R for a =1

1
R

 

 

 

 

 

 

 

 

 

 

 

 

 

(e) J/st vs. R for 0 =1
1.05 . . . . 5 :
1 ‘ 5'
.1 . ,1
$0.95 5
“E. .' —a = 0
5 ---a = 1
0.9; .f
0'851 1.2 1Z4_1I6 1:8 2 00 2
90%)
(h) T* vs- 10 = 90%)

(g) Pf/pgat vs. R for a = 1

Figure 5.5: Panels (a) - (h) correspond to figures 1 - 8 of [48] provided that the

boundary condition 7'0 = 5R0 is replaced by the boundary condition orr(ro) = 0.
This is the problem considered in Section 5.3. Here the effect of increasing the nondi-
mensionalized moment T* is to drive out the previously imbibed fluid. The ratio

p/ pg at is a formal calculation based on (5.22).

95

 

arr stress distribution.

5.4 Twist Induced Fluid Absorption

We now consider what can be viewed as the converse to the problem considered
in the previous section. That is r0 is now held at the free swelling value and the
twist applied at r = r,- such that the inner surface is free of normal traction. Thus

boundary conditions (5.19) are to be replaced by
Orr-(Ti) = 0, r(Ro) = CRO. (5.24)

Results associated with the solution to this problem are shown in Figure 5.4. In
particular, panels (a) - (g) of Figure 5.4 show the effect of the modified boundary
conditions (5.24) as they compare both to the case with boundary conditions (5.19)
as shown in Figures 1-7 of [48], and to the case with boundary conditions (5.23) as
shown in panels (a)-(g) of Figure 5.3.

Under boundary conditions (5.24) one finds that increasing 8 tends to cause
additional fluid uptake as shown by the decrease in r,- = r(Ri) in Panels (b) and (f)
of Figure 5.4. As in the previous cases, twist also causes fluid to redistribute from
the inner to the outer portion of the tube. The reader may have observed that we do
not yet give the torque vs. twist relation, i.e. the equivalent of Figure 8 of [48], and
hence the equivalent of panel (h) of Figure 5.3. A detailed discussion of the torque
vs. twist relation will be given shortly.

Figure 5.8(a) and 5.8(c) show the increase in V/st as 11) increases. Also Figures
5.8(b) and 5.8(d) give the arr stress distribution. It is to be noted that on- is no
longer monotonic in R for sufficiently large T*. In Figure 5.4 and Figure 5.8, for the
case a = 1, we use the same values T* = 1, 2,3,4 as in Figure 5.3, but we do not use

T* = 5. For a = 0 we must employ the very much lower values T* = 0.2, 0.4, 0.6, 0.8.

96

 

 

 

1.5 1 5

 

V/st

0.5*

----------------------------

 

 

 

 

 

 

 

 

 

 

0 1‘0 20
1/1

(a) V/st vs. 1,!) in radians for a = 0 (b) 0,»,— vs. R for a = 0

1.2 ' ' 1

 

 

1'

_o
oo

 

V/ View
9
O)

 

 

 

 

 

 

 

 

 

 

(c) V/st vs. 2,!) in radians for a = 1 (d) arr vs. R for a = 1
Figure 5.6: Additional graphs for the boundary value problem described in Figure

5.3. The dashed lines in panels (a) and (c) give the baseline associated with a pure
polymer (no liquid) by showing the volume ratio Vpoly/st = Jfll.

97

 

 

 

--------------
..—

........
‘‘‘‘‘‘‘‘‘‘‘‘‘
""""""
""""""""""
......
.....

 

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. .f
: ----- T" =OZ60 09" .+
0 9 . + 7“ = 0.80 g" ‘
' 1 1;5 2 1 1;5 2
R R
(a) J/st vs. Rforazo (b) f/R vs. Rfora=0
1.04
1.03'
1.02’
“as:
$1.01'
0...,
Q. _-
1.. .-
0.99"
0.98 * 1 1:5 2
1 1R" 2 E
(c)pf/p£atvs.I-2fora=0 (d)g=6——8vs.Rfora=0

98

 

 

 

1.3 . 1.05 +
, T“ = 0.000
---T* =1.000

 

1.2' ..... Tl ___ 3.000

 

 

 

J/Jf,
r/(CR)

......................
.......
‘
'0
‘I
I

..............................................
u
.......

 

 

 

 

 

 

 

0.9

 

(e) J/st versus R for (1 =1
1.15

 

 

 

 

 

 

 

R
(g) pf/pgat vs. R for (1 =1

Figure 5.7: Panels (a) - (g) correspond to figures 1 - 7 of [48] provided that the
boundary condition 1',- = (R,- is replaced by the boundary condition 0'1","(7'i) = O.
This is the problem considered in Section 5.4. The effect of increasing T * is now to

cause more liquid to enter into the system. The ratio pf / p}: at is a formal calculation
based on (5.22).

99

The reason for this will be discussed shortly. In addition, since fluid must pass from

the exterior liquid into the gel mixture, there is the possibility that the mixture could

become nonsaturated if there is insuflicient fluid in the bath to provide for the swelling

volume increase.

1.2

 

1.15’

1.1

V/ st

1.05- *

 

 

 

10 05 1
11)
(a) V/st vs. 21) in radians for a = 0
1.2 -

 

1.15*

1.1-

V/st

1.05-

 

 

 

10' 05 i
1P

(c) V/st vs. 112 in radians for a = 1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.3 .
—T* = 0.00
"'T" = 0.20
_ -------"T = 0.40 ,,,,
0.2 /- m-‘T = 0.60
g + T‘ = 0.80
i f
g: 0 1 f """"""""""""""""""""""" ~
b o: iii
0'52"; -------------------------
-0.1 -
1 1._5 2
R
(b) an- vs. R for a = 0
2.5
2‘ 1"“ -,
/ —T" = 0.00 W
3 ---T* = 1.00
1.5’ f ....... T! = 2.00 '
=1 f ----- T‘ = 3.00
\E: 1. z + T‘ = 4.00
b : ________________ .-.-. ...........
0.5*: .’.Iuln I
h’----_ ----------------
o . ______
-0.5 ‘
1 145 2

(d) arr vs R for 0 =1

Figure 5.8: Additional graphs for the boundary value problem described in Figure

5.4.

100

5.4.1 A Prescribed Inward Displacement of the Inner Radius

Before considering the possibility of a loss of saturation, we return to the issue of
why the maximum value of T* considered in Figure 5.8(b) is T* = 0.8 for a = 0 and,
similarly, why the maximum value of T* in Figure 5.8(d) is T* = 4.0 for a = 1. Let
us first consider the case a = 0. In this case our numerical method is unable to obtain
solutions if the twisting moment T" is too large. To explore this phenomenon consider
yet another boundary value problem where the boundary condition arr(r,-) = 0 is
replaced by a condition of prescribed displacement at the inner radius. Thus boundary

condition (5.19) is now replaced by

r02.) = 6412., r02.) = 4a., (5.25)

where 6 = Ti/ (CRi) is prescribed. Thus boundary condition (5.19) is the special case
for which 6 = 1. Our interest is in the case 0 < 6 < 1 so that the inner radius
displaces inward meaning that fluid is drawn into the cylinder.

It is also convenient to specify the torque T* in place of the twisting angle 11).
Figure 5.9(a) shows graphs of arr(r,-) vs. 6 at various fixed values of torque T*. For
6 = 1 we recover the result that T* 7E 0 implies arr(r,-) < 0. It is found that the
value of 01-1- at the inner radius is increasing as 6 decreasas from 6 = 1. A solution for
the problem characterized by the previous boundary condition (5.24) corresponds to
the locations where the curve in Figure 5.9(a) crosses the line orr(r,-) = 0. We find
that such a crossing takes place only if T* is sufficiently small, namely if T* < 0.937.
However it is found that if T* > 0.937 then such a crossing does not occur because
the curve of arr(r,-) vs. 6 is first increasing and then decreasing (as 6 decreases from
one) in a way that keeps the maximum value of 0'7")"(Ti) < 0. This is shown in Figure
5.9(a) where the T* = 1 curve remains below the arr(r,-) = 0 axis. Thus T* = 0.937

is a critical value of twisting moment beyond which we do not satisfy the condition

101

orr(r,-) = 0 for any 6 and hence do not obtain solutions for the previous boundary
conditions (5.24). Note from Figure 5.9(c) that if T* = 0.937 then the graph of
arr(r,-) vs. 6 will have a single tangential intersection with the line (rm-(ri) = 0.
Furthermore the graphs in Figure 5.9(c) show that if T* is slightly less than the
critical value 0.937 then the curves in Figure 5.9 intersect the Urr(7‘i) = 0 axis in two
locations, indicating the existence of two solutions to the problem with the previous
boundary condition (5.24). Here we remark that the graphs displayed in Figure 5.4
for T* = 0.2, 0.4, 0.6, 0.8 are associated with the intersection that is closest to 6 = 1.
For example, the T* = 0.8 graphs in Figure 5.4 (a) - (d) are associated with point A
on Fig 5.9(c). We shall refer to these as standard solutions for boundary conditions
(5.24).

A similar situation obtains when 0: > 0 in that there is still a critical value
of the twisting moment at which the graph of arr(r,-) vs. 6 has a single tangential
intersection the horizontal line arm-(ri) = 0. The critical twisting moment T* varies
with a. For 01 = 1 the critical twisting moment is T* = 4.041, which is why we

cannot take T* = 5 in Figure 5.4 when a = 1.

5.4.2 Further Consideration of the Problem with Boundary
Condition (5.24)

It follows that, in general, solutions to the problem with boundary conditions
(5.24) do not exist if T* exceeds its critical value. Conversely, if T* is less than its
critical value, then more than one solution may exist for the given T*. In addition,
for any solution to the problem with boundary conditions (5.24), one may calculate
6 2 72/ (CR1) and so identify the solution with a solution that employs boundary
conditions (5.25). In this way we have identified the standard solutions discussed in
Section 3.3.1 with special values of 6. Namely, these are the values of 6 for which

the intersection of the line 07.7.(ri) = 0 with a graph from Figure 5.9 occurs closer to

102

'r—__.w

 

 

6x 106dyne/cm2 , x 1,08 dime/077?2

 

 

 

 

 

 

   

 

 

 

 

 

 

2
—T* = 0.00
---T‘ = 0.20
4 ------- T" 2 0.40 l
~-~*‘T = 0.60
_ “T = 0.80
4: 2* ,,,,,,,,,, T* = 1.00 ‘
3:, ....
6
-2 -6
_4 , , 8 . . . 4;
0.6 0.8 1 0.2 0.4 0.6 0.8
6 6
(a) (b)

2

 

x 106dyrw cm

arr (Ti)

 

 

 

Figure 5.9: arr(r,-) versus 6 for or = 0. Panel (a) is the magnified portion of panel (b)
near 6 = 1. Panel (0) provides additional curves near 6 = l with the top curve for
T* = 08; bottom curve for T* = 1.0, and with additional curves for T* increments
of 0.01 in between.

103

6 = 1 than any other intersection.

We now turn our attention to the other values of 6 at which (rm-(n) = 0 when
there are multiple such intersections. For example, location B in panel (c) of Figure
5.9 corresponds to this solution for the case T* = 0.8 when a = 0. This also pro—
vides a solution to the boundary value problem with boundary conditions (5.24) and
corresponds to a smaller value of 6 and a larger value of «I; then that of the stan-
dard solution. We shall call these post-maximum solutions for reasons that will soon
be apparent. These post-maximum solutions are increasingly difficult to calculate
numerically as we encounter certain numerical instabilities as 6 gets small.

For a = 0 it is found that the post-maximum solutions can be calculated rela-
tively easily for 6 > 0.3. However for 6 < 0.3 they are highly sensitive to the initial
guess used in the iteration procedure. Furthermore for 6 less than about 0.15 we
cannot easily find an initial guess providing a converging solution. Indeed the min—
imum value of 6 for which we can find a converging solution is also dependent on
the value of T*. The curve of T* vs. ’l/J is shown in Figure 5.10(a). The absence of
solutions for T* greater than the critical value causes the graph of T* vs. a to exhibit
a local maximum at the critical value of T*. The end of the curve in Figure 5.10(a)
at «p = 1.7 represents the limit of computable solutions by our numerical method.
The ratio V/st is displayed in Figure 5.10(b) with the same range of 1/2 as in Figure
5.10(a).

A similar maximum twisting moment condition has been implicated in boundary
value problems in nonlinear elasticity for azimuthal shearing deformations of a form
similar to (5.5), but without a notion of swelling, when subject to boundary conditions
of fixed radial displacement in place of the boundary conditions considered here.
In [81] a phase plane analysis of the governing ordinary differential equations for
such a boundary value problem, using a compressible hyperelastic Blatz-Ko material

model, is found to generate a limiting moment condition for the construction of certain

104

 

 

 

 

 

 

 

 

1 - 1.4
l ordinary E post-maximum
0.8’ 1 solutions : solution
3 1.3' *— g__, l
0.6- 5 ,,, §
*8. ‘_ 5___> e12 5
M a . > :
ordinary : post-maximum :
solutions: solutions 1 1. i
0.2- : ~ ° 5
5 ¢=QW : w=0m
0 ‘ = .
0 1¢ 2 1o 1,, 2
(a) (b)

Figure 5.10: (a) Torque vs. twist relation with oz = 0 for the boundary conditions
(5.24). Here T* = 0.937, 11) = 0.90 corresponds to the local maximum. Also T* =
0.377, 1,6 = 1.70 is the effective limit of our numerical calculation and corresponds to
6 = 0.15. (b) Volume ratio V/st versus twist angle 1,0 for the same problem.

smooth deformation solutions. In [82] the numerical treatment of such a boundary
value problem for a broader class of Blatz-Ko models gives rise to numerical difficulties
if the applied twisting moment becomes too large, which is also suggestive of an
underlying maximum moment condition.

Figure 5.11 shows the numerical results corresponding to these post—maximum
solutions. The behavior of the curves in Figure 5.11 is similar to those shown in
panels (a)-(d) of Figure 5.4. Thus, at least with respect to these graphs, the post-
maximum solutions are not strongly distinguished from the standard solutions. Fluid
redistributes from the inner portion of the tube to the outer portion of the tube as T*
decreases. This, however, still corresponds to increasing a. In view of the decreasing

graph of torque vs. twist for the post-maximum solutions, one would also expect basic

105

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.5 1 _
0 _ 8 . x x': ,
S 1 $0.8"
a —T* = 0.600 t
---T* = 0.700
------ T* = 0.800 0.4,-1,3
----- T* = 0.900
~ T* = 0.937
0.01 145

  
   

R
(a) J/st vs. R for a = 0
1.1

 

 

 

 

 

 

 

 

1.
«$9:
$0.9
1...
Q
0.8
0.71 145
(c)pf/p£atvs.Rfora=0 (d)g=9—6vs.Rfora=0

Figure 5.11: The post—maximum solutions for a = 0. Thace correspond to a. contin-
uation of panels (a) - (d) of Figure 5.4. Now T* decreases with increasing 112 as is
evident from panel ((1).

differences with respect to any consideration of energetic stability (viz., eg., [83]).
5.5 Fluid Uptake Leading to Loss of Saturation

In the boundary value problem with condition (5.24) we have found that the
overall volume increases with increasing 1/2. All such volume increase is associated
with uptake of liquid, and it has so far been assumed that there is always enough
liquid available to support this volume increase. In other words, the analysis of
Section 3.3 is predicated on the assumption that there is sufficient liquid to keep the
system in a saturated state. We now consider the situation in which the amount of
liquid is limited, and in particular assume that all the liquid is imbibed at some point
in the loading procass. The instant at which this occurs corresponds to a transition
from a saturated state to a nonsaturated state. Under further increase in 1/2 the

system remains nonsaturated and the overall cylinder volume remains fixed at the

value V = Vpoly + Vliq-
Since V = 7rCL(rg - r?) and Vpoly = 1rL(Rg — R3), the constraint (2.20) can be
written

«(L03 — r3) = «L023 — R?) + mg. (5.26)

In addition, (5.24) gives r0 = CR0 whereupon it follows that the constraint is equiv-

alent to a requirement that

 

1 W t
r,- = \/(2Rg — 5023 — R22) — fl 2 Tz(ncmsa ). (5.27)

There will then also be a loss-of-saturation value of twist 1,!)(1’03) and a loss-of-
saturation value of moment T*(LOS).

In view of the additional specification of r,- in (5.27) the boundary conditions

107

 

 

(5.24) are, after loss of saturation, replaced by
arm-(6) = 0. r(R.) = rlnmatl, rate) = (BO. (528)

In contrast to the boundary conditions for the saturated cylinder, in which (5.24)
specified exactly one condition at both the inner and outer radii, now two conditions
are specified at the inner radius while one condition remains at the outer radius. It
is however possible to meet all of the conditions (5.28) because the uniform pressure
p in (2.22) is now added to all the normal stresses in (5.7) - (5.9). Note that this
constant term is differentiated out of the equilibrium equations and so (5.11) and
(5.12) remain as governing equations. These governing equations can now be solved
subject to (5.20) and to the second two conditions in (5.28). This would seem to be a
well posed problem, at least in terms of the number of boundary conditions that are
employed. Indeed it is similar to the problem associated with boundary conditions
(5.19) and to the problem associated with boundary conditions (5.25). After such a
solution is obtained, the value of p may then be selected so as to satisfy the remaining
condition orr(r,-) = 0 in (5.28).

As an example take a = 0 and suppose that the transition from a saturated
condition to a nonsaturated condition takes place when T* = 0.6. There is both a
standard saturated solution and a post-maximum saturated solution associated with
T* = 0.6. The former is with 1/) = 0.418 and 6 = 0.946; the latter is with 1,6 = 1.412
and 6 = 0.302. These correspond to different amounts of available liquid. Specifically,
loss of saturation occurs for the regular solution at T* = 0.6 if Vlz‘q = 40.469LRg. Loss
of saturation occurs for the post-maximum solution at T * = 0.6 if Vqu = 53.371LR‘E.
Alternatively stated, let us consider R0 = 2R,- and a = 0 for the two cases Vliq =
40.469LRS’ and Vliq = 53.371LRg. Then the boundary value problem with boundary

conditions (5.20) and (5.24) causes complete fluid absorption under increasing :6 for

108

(nonsat)
i

the case Vqu = 40.469LR? when 21) = 0.418, in which case r = 0.946(Ri; and

causes complete fluid absorption for the case Vlz'q = 53.371LR? when w = 1.412, in

(nonsat)
i

which case r = 0.302CR,-. The two values of Viz-q are here specially chosen to
give the same value T* = T*(LOS) = 0.6. Thus the two cases correspond respectively
to a standard saturated solution transitioning to a nonsaturated solution and to a
post-maximum saturated solution transitioning to a nonsaturated solution, both for
the same value T*.

The boundary value problem for each of the two cases with T* (L05) = 0.6 is
solved numerically by the methodology outlined above. Since the boundary condition
arr (7‘2“) = 0 is met by the freedom to choose 19, the solution method is the same as
that in which the radial displacement is specified at both R,- and R0.

Figure 5.12 shows the mechanical behavior of the nonsaturated states associated
with the transition taking place at the standard solution discussed above, i.e. qu =
40.469LR9, T* = 0.6, w = 0.418. Panels (a) - (d) of Figure 5.12 should be compared
with panels (a) - (d) of Figure 5.4 since Figure 5.4 shows the mechanical behavior
for the same system prior to the loss of saturation. In particular the graphs for
the transition value T* = 0.6 are the same for both figures. In the context of this
example, the T* = 0.6 curve in Figure 5.4 should now be disregarded since there is
now insufficient fluid for a saturated solution when T* = 0.6. In a similar fashion,
panel (e) of Figure 5.12, which gives arr(r), should be compared to panel (b) of
Figure 5.8. As T* increases one again finds that fluid redistributes from the inner
portion to the outer portion (panels (a) and (c)). Unlike the saturated solutions prior
to transition, now both r,- and r0 remain fixed as T* increases (panel (b)). Also, the
twist angle 16 continues to increase with T* (panel (d)).

Figure 5.13 show the mechanical behavior of the nonsaturated states associated

with the transition taking place at the post-maximum solution discussed above, i.e.

V“, = 53.371LR§, T* = 0.6, w = 1.412, 6 = 0.302. Hence the T* = 0.6 curves in

109

 

 

0.98» ]

 

“00.96

 

 

 

 

 

 

 

 

 

 

 

 

0.94-“""'
1 1; 2 0'9‘1 L5 2
R R
(a) J/st vs. R (b) r/R vs. R
1.5 - 2 =x106 dyne/cmz
.0 '

 

 

 

 

 

 

 

0 - 5
1 1.5 2 ‘
R" 01 1.71.5 2
(c)g=0—6vs. R (d)01-rvs. R

Figure 5.12: Mechanical behavior for the nonsaturated solutions when loss of sat-
uration takes place at a standard solution with T* = 0.6. Here a = 0 and
Vlz'q = 40.469LRg. The solid curve in each panel corresponds to the instant of
transition while the other curves correspond to increasing T*. The values of 1,!) asso-
ciated with T* 2 0.60, 0.80, 1.00, 1.20, 1.40 are 111 20.418, 0.562, 0.708, 0.858, 1.011
respectively.

110

 

panels (a) - (d) are identical to those displayed in panels (a) - (d) of Figure 5.11. The
curves in Figure 5.11 give the behavior prior to the loss of saturation. The overall
behavior displayed by the curves in Figure 5.13 are similar to those in Figure 5.12.
In particular 1,!) increases with T * after loss of saturation. Thus the negative slope
of the torque-twist graph for the post-maximum solutions prior to loss of saturation
immediately changes to a positive slope torque-twist response after loss of saturation.

Figure 5.14 shows the nonsaturated torque-twist response curve corresponding to
each of the two transitions for T* = 0.6. These curves are connected to the previously
displayed response curve of Figure 5.10. In each case, loss of saturation is associated
with an abrupt stiffening of this macroscopic response. In particular, when loss of
saturation occurs at the post-maximum solution, the torque - twist rasponse is once
again monotonically increasing.

As an additional comparison consider again the same standard solution transition
state as that in Figure 5.12, namely the case with Vh-q = 40.469LRS’, so that transition
occurs at T * = 0.6, 21) = 0.418 with 6 = 0.946. Specially, for this Vh-q = 40.469LR;3
consider the three twist values 11) = 0.220, 0.418, 0.620. Figure 5.15 comparas
the associated sequence of mechanical response curves for the saturated, transition
and nonsaturated states that correspond respectively to 1,6 = 022,212 = 0418,11) =
0.62. In addition, this figure also shows the saturated response when 11; = 0.62
(which can only occur for a larger value of Vlz‘q)- A comparison of the saturated and
nonsaturated response for 2/1 = 0.62 makes evident the continuous volume increase
and inner diameter reduction for the saturated case. In contrast, once the system
is nonsaturated the value of Ti remains fixed (panel (b)) and the overall volume is
unchanged. In particular, the curves in panel (a) for w = 0.418 and for 1,0 = 0.62 after
loss of saturation cross over each other so as to permit J to integrate up to the same
value. Also, for r > r,, the stress Urr is everywhere greater for the nonsaturated case

with 1,0 = 0.62 than it is for the saturated case with this same 1/1.

111

 

 

 

   

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 1
0.8
. E0 6 .1333}?!
—T* = 0.60 E ' —T: = 0-60
---T* = 0.80 "'T = 0-80
------- T* = 1.00 . 0 4» ,s.-’.-;I>" T * = 1-00
..... T* = 120 I ’ ""'T* = 1.20
. T* = 1.40 * T = 1.40
01 L5 2 0'21 15 2
R R
(a) J/st vs. R (b) r/R vs. R
5 . 8X107dyne/crn2
61
$4“
------------------------- l
0 - s”
1 5 2 0 ‘
1 1.5 R 2

 

 

1 I

R _ _

(c)g=0—evs. R (d) arr vs. R

Figure 5.13: Mechanical behavior for nonsaturated solutions when loss of saturation
takes place at the post-maximum solution with T* = 0.6. Here a = 0 and qu =
53.371LRE’. The solid curve in each panel corresponds to the instant of transition

while the other curves correspond to increasing T*. The values of 1]) associated with

T” = 0.60, 0.80, 1.00, 1.20, 1.40 are 111 = 1.412, 2.053, 2.800, 3.632, 4.533 respectively.

112

 

T*

0.5-

 

 

 

%061163258354

Figure 5.14: Saturated and nonsaturated torque—twist response for a = 0. Two
examples are shown, each corresponding to T* = 0.6. The square markers on the
nonsaturated curves represent computational points and the intervening dots are to
guide the eye. The transition on the left occurs at a standard solution (corresponding
to the solutions displayed in Figure 5.12). The transition on the right occurs at a
post-maximum solution (corresponding to the solutions displayed in Figure 5.13).

113

 

1.15

 

 

 

‘n‘ I ‘ . ' ‘3‘.
.“fl '.‘|-‘ ” ""~ .3...
...... 0.95 pf *
V. 0'. Q.
a) 1 .05 ’ , """""""""" ‘ A
“ u ' ' - :-
«I, "1' Q:
3 .3 ”:" V
\ J. ’ f v
- I 6
l 5 3’ \
r I k
1 K
I : .l 0 9 .
I 5 I I

 

_i —16 = 0.22(sat)
0.3595 ‘1‘, ---16 = 0.42(tran) '

 

 

 

 

 

 

 

 

 

 

 

:3" ------- 16 = 0.62(sat)
O 9,’ ----- 161: 0.62(non) 0 85 .
' 1 L5 2 ' 1 L5 2
R _ _ R_
(a) J/st vs. R (b) f/R vs. R
5
8X 10 dyne/cm2
"""""""""" l

 

 

 

 

 

 

 

 

R 2
(c)g=0—8vs. R (d) an- vs. R

Figure 5.15: Saturated, transition and nonsaturated states associated with loss of
saturation taking place from a standard solution when 1,6 = 0.42. Thus 16 = 0.62
corresponds to a nonsaturated state. For comparison, the response for a hypothetical
16 = 0.62 saturated state is also shown.

114

In a similar fashion, Figure 5.16 provides the same type of comparison as in Figure
5.15 for the case in which loss of saturation occurrs at a post-maximum saturated
state. Here it is to be noted for 16 = 1.6 that the saturated case is more depleted of
fluid near r = r,- than is the nonsaturated case. In addition the a" profile for the
16 = 1.6 nonsaturated case is no longer everywhere above the 16 = 1.6 profile for the
saturated case. The exceptional regions where this ordering in arr ceases to hold is

again confined to a zone near r = 7’1-

5.6 Summary

In this chapter, the problem of twisting a swollen tube is formulated based on
the theoretical model described in Chapter 2. We started by replicating Wineman
and Rajagopal’s work [48] where they consider the problem of a volume preserved
swollen cylinder under twisting in a saturated state. We obtained the exact same
graphs as in [48]. Then more generally we consider different boundary conditions that
do not necessarily conserve the total amount of liquid in the cylinder, we discover
additional interesting effects when there is relative twist 16 7t 0 between the inner
and outer radii. Replacing the condition that the outer radius maintains its free
swelling value by a condition of zero normal traction causes this outer radius to
displace inward, thus decreasing the overall volume. This requires fluid to exit the
system. Conversely, replacing the condition that the inner radius maintains its free
swelling value by a condition of zero normal traction causes this inner radius to
displace inward, thus increasing the overall volume. This requires fluid to enter the
system. In all cases we find for 16 aé 0 that the fluid concertration in the cylinder is an
increasing function of the radius. As was the case in [48], the amount of redistribution
increases with I16] for the broader class of problems considered here. In addition, for

the problem characterized by the zero normal traction condition at the inner radius,

115

 

 

    

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

9",”
1 5 fi . 0.8- 4.6”
693“" 1,60?
.7; 65: 0.6' 5
’5 1 :1; 25-"
E F 5:":
'7: —16 = 1.21(sat) 0'4/1, —16 = 1.21(sat)
0 56' "'1/1 = 1410M”) , :5“~ ---16 = 1.41(tran)
' ------- 1) = 1.60(sat) 0.25 ....... ¢ = 1.60(sat)
""" 1b = 1150010”) .-.-.¢ = 1.60(non)
01 1_:_5 2 01 145 2
R _ _ R_
(a) J/st vs. R (b) f/R vs. R
2 8X 106dyne/cm2
1.5 "A...“.;'.".‘:'.‘:‘.".-.-.-."-'-'-'-“"‘ “m“ ‘

 

—16 = 1.21(sat)
---16 = 1.41(tran) .
------- 16 = 1.60(sat)
----- 16 = 1.60(non)

1:5 2
R
(c)g=0—9vs. R (d) 07-,- vs. R

 

 

 

 

 

 

 

 

 

 

 

Figure 5.16: Saturated, transition and nonsaturated states associated with loss of
saturation taking place from a post-maximum solution when 16 = 1.41. Thus 16 = 1.60
corresponds to a nonsaturated state. For comparison, the response for a hypothetical
1,6 = 1.60 saturated state is also shown.

116

we find that the relation between the twist angle 16 and the twisting moment T may
lose monotonicity so that there is a maximum twisting moment that can be sustained
under such circumstances.

Furthermore, for those situations that give rise to fluid uptake, we consider the
consequences if, at some point in the process, all of the fluid is drawn into the system
so that no additional uptake is possible. Then the transition from the saturated state
to nonsaturated state occurs. The overall volume of the cylinder cannot change once
the system loses saturation. The result shows that the now fixed amount of fluid
continues to redistribute as 16 changes. It is also found that the relation between the

twist angle 16 and the twisting moment T abruptly stiffens after loss of saturation.

117

 

Chapter 6

Seepage Through a Fiber
Reinforced Hyperelastic Layer

The effect of a fluid solvent on the swelling of a rubbery polymer is of longstand-
ing interest [44]. In particular, extensive experiments have been performed on the
pressure-induced diffusion of fluid through polymer networks [20—23, 84,85]. One of
the heavily referred set of pressure-induced diffusion experiments was conducted by
Paul and Ebra-Lima [24]. They recorded the nonlinear relation between the fluid
flux and the pressure difference driving the solvent through a highly swollen poly-
mer membrane layer. They also presented a combined thermodynamic and diffusion
theory for this pressure—induced transport. In addition, uniaxial stretch experiments
were made on the gum rubber material immersed in twelve different organic solvents
that were tested in the diffusion experiments. The elastic and thermodynamic con-
stitutive parameters were determined from these tests by using the Flory-Huggins
theory. These parameters were then used widely in later works that seek to model
such pressure-induced diffusive transport [46—49]. It is to be remarked that although
we use the term ”diffusion transport” here, other intellectual traditions may prefer
the term ”seepage”. Here we use there terms interchangeably.

A broader continuum mechanical framework can be invoked to deal with sepa-
rate mechanical balance principles for both solid and fluid component. This broader

framework is known by a number of names, including: large deformation mixture

118

theory, the theory of interacting continua, and the large deformation biphasic theory.
Examples of mixture theory type treatments applied to problems involving these sys-
tems include [47 , 71, 86], where a major focus is on time dependent swelling as liquid
diffuses within the elastomeric matrix.

In this chapter we consider similar seepage problems for both isotropic and
anisotropic (fiber reinforced) gels. In particular we shall use a framework recently
considered by Pence [87] that provides for consistency with equilibrium problems in-
volving solvent redistribution such as considered earlier in this thesis and also by

other researchers.

6.1 A Biphasic Mixture Theory

6.1.1 Kinematics and Notation

Following Shi et al. [46], we let Xf and X3 denote reference positions for fluid
and particles respectively in the reference state prior to mixing. Positions of these
particles at time t are determined by the respective functions xf (Xf , t) and x3 (X3, t)
(Figure 6.1). We then let x denote a location in the current configuration at time
t, it is assumed that each such spatial location is occupied by both fluid and solid

particles.

The fluid and solid velocity field are then given by

f a
f=__aX s=__X 61
v at xf’ v 61 X3' (')

The deformation gradient tensor associated with the motion of the solid constituent

is given by

3

(9x

F=aX3.

 

(6.2)

119

Solid

Fluid Fluid Particle

Particle

   
  

    
     

Particle

‘x3

Solid

Particle
X 2 X x2
1
X1 Reference Current
Configuration Configuration

Figure 6.1: Mapping of fluid and solid particles

J = detF, (6.3)

then J is the ratio of the current volume occupied by a collection of solid particles to

the volume that these same particles occupy in the reference configuration. Let pf

denote the fluid density in the overall mixture, i.e., the mass of the fluid constituent

per unit mixture volume in the current configuration. Conversely, let pg denote the
f

fluid density in the reference state. This p0 is a constant, whereas pf is a function of

x and t. The fluid density is then subject to the mass conservation equation:
iap—f+V-(pfvf)=0. (6.4)
at .

The solid densities p3 in the current configuration and pa in the reference configuration
are defined in a similar fashion. Then the mass conservation of the solid phase is
expressed as

I
s_ .9
p — JPO' (6.5)

The assumption that the volume of the mixture is the sum of volume of fluid

120

and solid constituents prior to mixing gives

_ . (6.6)
P0

By substituting from (6.5) into (6.6) one obtains

pf = (1 — :1?) pg. (6.7)

Since the total mass is the sum of the fluid and solid mass, the overall density p is
given by

1 1
p=pf+ps=7p8+(1-j)pg-

(6.8)
The seepage velocity w is the velocity of the fluid with respect to the solid, i.e.,

w = vf — vs. (6.9)
The final kinematic quantities to be introduced are the velocity gradient tensors
associated with the respective motion of the fluid and solid constituents, namely,

Lf=V®vf, L3=V®v3.

(6.10)

6.1.2 Fluid and Solid Partial Stresses Under Volume Con-
straint

In the absence of external body forces, and neglecting inertial effects, the equi-

librium equations of the solid and fluid phase are given by

V~03+b=0, Voaf—br—O,

(6.11)

121

 

where as, of and b are, respectively, the solid partial stress, the fluid partial stress
and the possible dissipative interactive force between the solid and fluid phase. The

tot

total stress a = 0'3 + of therefore obeys

v - at“ = 0. (6.12)

Let \Ilcur = \Ilcur(F) be the Helmholtz free energy of the mixture per unit mass in
the current configuration. Based on a mixture theory treatment involving the first
and second thermodynamic laws (details can be found in [87]), the fluid and solid
partial stresses and the interactive force with the constraint of volume preservation

(6.6) are given by

 

b = pfaq'm :VF — p—lvpf + 13 (6.13)
6F f
P0
as = page—EFT — 111 + 63, of = —£f-p11 + (if. (6.14)
6F J 25

Equation (6.11) are standard in the theory of interacting continua [48] whereas equa-
tions (6.13) and (6.14) generalize the isotropic theory equations (9), (10) and (11)
of [48]. Here 191 is a Lagrange multiplier which arises due to the volume constraint
(6.6). The tensors 63 and &f are possible dissipative stresses and f) is a possible

dissipative interactive force. These are required to obey
asstzo, aszfzo, B-wzo. (6.15)

It follows from (6.11)2, (6.13) and (6.14)2 that

f
Vpl = waved... + £9 (v . 6f + 13) . (6.16)

122

Eliminating p1 in (6.13) and (6.14) with the use of (6.7) and (6.11) gives
(J —1)v- (amix + {73) = v - 6f — .113, (6.17)
where 0'77”.“c is defined by

' \I/

Exactly the same as in (2.11), let W be the Helmholtz free energy of the mixture
per unit volume of the solid in its reference state. Therefore W is related to \I/cur,

the Helmholtz energy per unit mixture mass in the current configuration, by

 

p3 1 1 3 W
W=PJ‘I’cur=P—g‘pcur fi’ ‘I’cur=—J'W=—p—S = f.
p ppo 125+ (J —1)p0
(6.19)
It then follows that (6.18) is alternatively expressed as
l 8W T
mm: _ _
_. J 8F F . (6.20)

If {73, &f and f) all vanish, we then obtain from (6.17) the equilibrium equation
V 10mix = 0. (6.21)

To also account for the dissipative affects we will use (6.17) with &f = 6'3 = 0 but
b aé 0 as the governing equation for the fluid seepage problems that will be introduced

in the next section.

6.1.3 Boundary Conditions

The boundary of the mixture can admit a variety of situations including: (a)

contact with the pure fluid phase, (b) contact with a rigid support structure that

123

 

strictly isolates the mixture from any pure fluid phase, and (c) contact with a porous
rigid support that permits interaction with a pure fluid phase on the other side of the
support. To specify appropriate boundary conditions, let n be the directional normal
vector to the mixture external surface, let Pf be the pressure of any pure fluid phase
that is adjacent to the mixture and let t be the traction provided by any solid support
structure. In particular, the traction t may or may not permit interaction with the

pure fluid. Certainly one boundary condition is

totn = ttot’

a' where ttOt = t — an. (6.22)

Note that t = 0 for situation (a) given above, whereas Pf = 0 for situation (b) given
above. For situation (c) given above both t and Pf would typically be nonzero. In
general (6.22) does not provide enough boundary conditions to solve all boundary
value problems of interest. Thus more boundary conditions are usually required.
There is some controversy as to what these additional boundary conditions should
be [88]. Here we briefly discuss two proposals for the additional boundary conditions,
one is called the traction splitting boundary conditions, the other is based on an
argument involving chemical potentials. After this brief discussion we will adopt the
chemical potential boundary condition in what follows.

The first proposal is the traction splitting condition. As originally stated in
[89] (see page 36-37), the traction splitting treatment states that the partial surface
tractions associated with the solid and fluid, respectively, are proportional to their

local volume fractions. In our notation this would be expressed as

f f 6V3 3 1
f ___ 5V ttot = p_ttot = _ 1 tot s = ttot : P_ttot = _ttot
” ” 6V pf 1 J t ’ a “ 6V p3 J ’
0
(6.23)

where n is again the unit normal vector. However later (see page 42 of [89]) it is

124

argued that the porous solid support structure has “little effect on the flow of the
fluid constituent” and this is used in [89] to restrict the splitting in (6.23)1 to only
the fluid pressure. Letting Pf be the fluid pressure this gives an addition boundary
condition in the form of = — (1 - §) Pf.

The second, and alternative, proposal for boundary conditions appropriate to a
mixture of the type considered here consisting of an incompressible solid and incom-
pressible fluid is often called the chemical potential continuity boundary condition.
This type of treatment can be found in [71] and [90] where the chemical potential
continuity boundary conditions are obtained through a variational approach of ex-
tremizing the overall free emery of the system. In the present notation, this boundary
condition is stated in the following equations (see (39), (40) and (22) in [71])

816 T ~
[EFF +(16-p6K)I[n=t—an, and K=Kf, (6.24)

where

l
Kf E 27 (Pf + #lz’q) . (6.25)
0

Here 16 is the system Helmholtz free energy per unit current volume, K = K (x) is
the chemical potential of the mixture, Pf is the fluid pressure at the boundary as
given earlier in connection with (6.22), K f is the chemical potential of the fluid, and
”fig is a constant that gives the Helmholtz free energy per unit current volume of the
pure fluid phase. It is to be clarified that 16 is defined as Helmholtz free energy per
unit current volume of the system, while \Ilcur in (6.13) and (6.14) is defined as the
free energy per unit current mass of the mixture, which does not take into account

the free energy of each phase ”polg) and ”liq‘ By this we mean

.. 1 l
w : p‘I’C’uT "l' (I. - j) “liq ‘1'“ jflpoly. (6.26)

125

 

Comparing the system free energy 16 with the system free energy E per unit

reference volume as defined in (5.14) it follows that

~ 1 1 1 1
16 = 71E = 7W + (1— j) #liq + j/‘poly- (6.27)

It then follows from (6.24) with the use of (6.27) that

l 3W T
(j—aFF + #liq) n = t + Mllqn, (6.28)

which immediately leads to

a ' n = t (6.29)
by the definition of am” in (6.20).
For our purposes in what follows we note that (6.17) with &f = 6’3 = 0 gives

V - ami115 + —b = 0. (6.30)

In addition, following Parasad and Rajagopal [88] the dissipative interactive force 6
will eventually be taken to be proportional to the seepage velocity w by virtue of the

following form
A a 1 n—l
= ._ __ . 1
b J (J _ 1) w, (6 3 )

where a > 0 is a material diffusive constant, and n > 0.

6.2 Free Swelling of a Fiber Reinforced Hypere-
lastic Layer

We consider a fiber reinforced gel consisting of an anisotropic hyperelastic solid

that is immersed in solvent. Prior to immersion the solid is taken to be in the form

126

of a layer that is reinforced with parallel straight fibers making an angle (2 with the
layer’ 8 thickness direction( (Figure 6. 2). The layer’s initial thickness is H. Let (X, Y,

ZZZZZZZZZZZZZZ/{Q Step1 ///////////////1‘}

Coordinate
Transformation
(Relabeling)

 

Free swelling that
preserves the
fiber direction

Step 2

 

 

Coordinate A
3; Transformation 9': 8
(Rotation) 7’)
it

Figure 6.2: Free swelling of a fiber reinforced layer

Z) denote the position of particles in the reference configuration. Let M denote the
unit vector of the fiber reinforcing direction in the reference configuration and assume

that the fibers in the reference configuration lie in the (X, Z)-plane. Thus we take

MzcosSleX+sinQeZ. (6.32)

We assume that the fibers are at their natural length in this reference configuration.
Now let the reinforced layer be immersed in a fluid environment and swell freely to
an equilibrated state with (it, 17, :2) describing position in the swollen configuration.

We wish to describe the resulting free swelling in a manner such that the (X, Y)-

127

 

plane has the same orientation as the (it, m-plane. To describe this homogeneous
deformation it is convenient to decompose the free swelling into three steps (Figure
6.2). The first step is a coordinate relabeling transformation that transforms the (X,
Y, Z) coordinates to (s, 17, c) by a rotation through angle (2 about the Y axis. This

aligns the 5 axis with the fiber direction. The mapping function and its gradient are

 

 

given by
€=XcosQ+ZsinQ cos!) 0 sinQ
n = Y , F1 = 0 1 0 . (step1)-
<=—XsinQ+Zcosf2 —sinfl 0 cosfl
’ (6.33)

The second step is a free swelling that preserves the coordinate directions (a, r), c).
Let A1 be the stretch ratio along the direction of the fibers while the stretches in all

the transverse directions are the same and so denoted by A2. This gives

 

 

83 = 2118 A1 0 0
6 = /\2€ 0 0 A2

The third step is a coordinate rotation transformation that rotates the coordinates

(6‘, 1), c“) to (is, 37, 2) by an as yet unspecified angle -w about the 7‘) axis

 

 

izécosw—é‘sinw cosw 0 —sinw
’9 = 0 2 F3 = 0 1 0 , (step 3).
2=ésinw+<fcosw sinw 0 cosw
’ (6.35)

Therefore the deformation gradient F associated with the mapping from the reference

configuration (X, Y, Z) to the free swelling configuration (:6, g), 2) is given by

128

A1 cochosw + A2 sianinw 0 A1 sichosw — A2costinw

 

 

F = F3F2F1 = 0 A2 0 ’
A1 costinw — Agsinflcosw 0 A1 sianinw + A2 cochosw
(6.36)
and so
i" = (A1 cochosw + A2 sianinw) X + (A1 sichosw — A2 costinw) Z,
= A2Y, (6.37)

2 = (A1 costinw — A2 sichosw)X + (A1 sianinw + A2cochosw) Z.

The left and right Cauchy-Green deformation tensor that follows from (6.36) is

 

 

given as
A2 cos2 w + A2 sin2 w 0 (A? — Ag) sinw cosw
B = FFT = 0 Ag 0 (6.38)
(A? — Ag) sinwcosw 0 A2 sin2 w + A2 cos2 w
and
A12cos2 Q + A% sin2 0 0 (A? — Ag) sin 9 cos 9
C = FTF = 0 Ag 0 . (6.39)
_(A112 — Ag) sin 0 cos 9 0 A? sin2 0 + Ag cos2 9—

 

 

Thus 11— _ tr(B)= A2 + 2A2, 12: % [tr(B)2 — tr(B2)] = 2A§A§ + X21 and J =
fl“ = detF = A93.

The angle w is now chosen so as to satisfy

A1 costinw — A2 sichosw = 0 z)» w = arctan (3:2 tan 9). (6.40)
1

129

With this choice of w the planes Z = constant remain at their original orientation

after free swelling. It then follows from (6.37) that

A Q A ' Q
1cos + 25m

 

 

 

 

 

 

 

x = A2 2 A2 2
1+X§tan0 1+chot Q
+ A1 8:19 _ A2 c359 Z,
A A
1 t 252 1 t2Q
\/ + A? an \/ + Sgco

3; : )‘2Ya (6.41)
2 = A1 st + A2 0080

 

 

 

 

A2 A2
1+—§cot2n 1+—§tan2a
"2 A1

The requirement (6.40) is simply one specification of the rigid body rotation
associated with a free swelling homogeneous deformation. For our purpose it is the
particular rigid body rotation that corresponds to the boundary value problem of
interest in this study. This is because it preserves the orientation of the planes
Z = O which will be a fixed boundary in what follows. Note also for the special case
A1 = A2 = A that (6.41) simplifies as expected to the equi-axial free swelling 3“: = AX,
3) = AY, 2 = AZ.

6.3 Constitutive Theory

The free energy function W will be taken as the sum of the three terms. These
are: an isotropic elastic energy <I>(Il,12) due to the deformation of the matrix in
which the fibers are embedded, the mixing energy H (J) due to the mixing entropy

and mixing enthalpy, and an extra term 7(1) fib(14) due to the stretching of the fiber

130

 

reinforcement. This gives
W = 9(11,12)+ H(J) + 7q>fib(14)- (55-42)

We take the Mooney-Rivlin model (2.13) for the elastic energy term <I> and the Flory-
Huggins form (2.5) for the mixing energy term H. For the last term of (6.42), 'y > 0 is
a material modulus and I4 2 FM - FM is the square of the fiber stretch as discussed
in [91]. On the basis of models in [91] and [92], the fiber energy <I> fib(14) in (6.42)

will be taken as

(172604) = $414 - 1)2- (6-43)

We first consider the homogeneous deformation for free swelling with F given by

(6.36). Then (6.21) is automatically satisfied and it follows from (6.29) with t = O

for all possible orientations n that 0min: must vanish identically. It thus follows from

(6.20) and (6.42) that

 

2 84> 2 6(1) 66» dH d<I>f.-b
——— —— — FM FM: 0. .44
JaIQB + J 2(61—1+ 11612) 3+ dJI+2J 414 ® (6 )

Expanding the above equations gives

,u (1 — g) (A2 cos2w + A3361? 6))
+115 [A3 31112 w + A2A§ (51624; + 2cos2 w)] + Jh(J) + 27 (A2 — 1) A1 cos2 w- _ o,
,u (1 —§+5A§) (A2 —A3) +27(A21’- 1) A2 = 0,
p [1 — 6 +5 (A22 + A2)] A2 + Jh.(J) = o, (6.45)
11(1—5) (A2sin2w + A2 2cos2 6))

+p§ [A2c032w+A¥Ag(coszw+2sin2w)] +Jh(J)+2'y(A2—1)A25in2w=0.

131

It can be shown that if (6.45)2 and (6.45)3 hold, then (6.45)1 and (6.45)4 are
automatically satisfied. The free swelling principal stretches A1 and A2 are now
determined from (6.45)2 and (6.45)3. Note that (6.45)2 and (6.45)3 are independent
of w. This is as expected since one may verify that (6.45)2 and (6.45) 3 follow from a
principal frame analysis of (6.44) by taking F = A21 + (A1 - A2) M (8) M.

The stretches A1, A2 and the corresponding values of J are plotted in Figure 6.3
for various values of '7. The special free swelling case for which '7 = 0 means that
there is no fiber effect so that A1 = A2 and J = A? Note if 'y > 0 then A2 > A1
since the additional resistance of the fibers causes less extension in the fiber direction
under free swelling. Note that the orientation of the fibers w after the free swelling as
given by (6.38) is dependent upon the reinforcement coefficient 7 since 7 determines
the difference between the fiber stretch A1 and the stretch A2 that is orthogonal to
the fibers. In particular, since A2 > A1, it follows from (6.40) that w 2 9 (w > 9
if 0 < 9 < 7r/2) and also that the difference between w and 9 increases with 7 (see

Figure 6.4).

6.4 Steady State Diffusion through a Fiber Rein-
forced Layer

Now we consider pressure—driven diffusion of fluid through such a swollen fiber
reinforced layer. The bottom of the layer is taken to be fixed to a rigid porous
plate that is attached to the gel after free swelling (Figure 6.5). Let (x, y, 2) denote
coordinates in the deformed configuration. The pressure driven diffusion is anticipated
to cause either a compaction or an expansion in the layer, and the presence of the fibers
is anticipated to cause a further shearing deformation. Based on (6.41) the overall

deformation from the original reference state to the deformed state is therefore taken

132

A1, A2 and J

A1, A2 and J

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

\
4-.. - "A1 .........................................
..... A2
—-J = A1A§
2 ~_£~_;;_~_-;_-_-;;-;-.-.~.-.-.-.-.~-.-.--.--.-.-~.-.-1.x; ;..._._.._._.._,..._._,..__...,_._._,..__
O
O 1 2
7/41
(a) g = o
4
3 ....... ---A1 ......................................
..... A2
2 ....... “J : 2‘12‘3 ......................................
1 ------------------
0 1 2
’Y/ 1a
(b) g = 0.5

133

 

 

 

 

 

 

'3 2.5- ................. _ .......................................................... .
73 ---A1

Cl

‘6 2- "2‘2 .........................................
,9 ——J = AlAg

2 1.5;;;.I,.._.k._._.‘.,_.,._.,_,,..,....-.-.. .... v.'.‘.-'.t"..'.'..'.'f’.‘.'.'..'.'". .'."'.'.'.'T'-..'.'.T'..'."

 

 

 

®)€=1

Figure 6.3: Stretch ratios of free swelling with fibers versus 7 for 5 = O, 5 = 0.5 and
5 = 1(M/p = 100, x = 0.425)

to be
:1: = :3: + 3(2) = A1X + f(Z), (6.46)
Z=Mfl=gwl (6%)

where
A1 = A1 cos 9 cosw + A2 sin Q sin a). (6.49)

Here the two expressions that multiply Z in (6.41)1 and (6.41)3 have been absorbed
into the functions f (Z) and g(Z) respectively. Note that Q is the original fiber
orientation and so known at the outset while A1 and A2 are regarded as known from
the free swelling analysis (e.g. from Figure 6.3). Thus a) is known in terms of this

parameter from (6.40) after which A1 follows from (6.49). Note also that certain

134

 

 

 

 

 

 

 

 

 

 
 
 
 

0.3'
(20.2- 5
l —9 = 0(1)
3 _ —Q = 7r/8 (2)
0'1 —o = ”/6 (3)
—Q = 7r/4 (4)
—Q = 7r/2 (6)
0 2 *1/14 4
(a) w—Q versus 7/11 with fixed!)
0.35
0.25r
C} 0.2”
| 0.15’
3
0.1“
0.05’

 

 

 

 

 

 

(b) w — Q versus (2 with fixed 'y/u

Figure 6.4: The variation of the change in orientation angle 11) — Q (in radians) with

é = 0 and M/y = 100, x = 0.425.

135

care must be taken with (6.46) since (6.37) shows that f (Z) 74 3(3) since, instead,
f(Z) = (A1 sichosw — A2 costinw)Z + 3(2).

 

 

 

 

   

Z Z = H
Y //
\
X Z = 0 fibers
(3.) Original reference configuration (prior to immersion in solvent)
_ PH 1
11:12:61211111
2 ivlvevvvwvivvv v freesurface
y 1
0)
A : __ fixed surface
x 1 ' 1 i 1 r l l g l l I l l 1

(b) Depiction of the praesure part of the boundary conditions with respect to
the free swelling configuration

Figure 6.5: Representation of a pressure differential that will cause diffusional fluid
seepage through a fiber reinforced layer

We first consider a steady-state diffusion in which case the solid particlas are not
in motion (v3 = 0). The fluid motion is assumed to be confined in the z direction,
meaning

vi = f(z)ez. (6.50)

Therefore the seepage velocity w = vf.
The deformation gradient tensor F and left Cauchy-Green deformation tensor B
associated with the mapping (6.46)—(6.48) from the reference configuration (X, Y, Z)

to the swollen configuration (x, y, 2) associated with the steady state diffusion are

136

then given by

A1 0 f’ A? + f’2 o f’g’
F = 0 A2 0 , B = 0 A? 0 , (6-51)
0 0 g, f/g/ 0 9/2

 

 

 

 

where f’ = df/dZ, g' = dg/dZ. Thus 11: A? +A§ + f'2 +g'2,12 = [fig + Agf’2 +
(A? + A?) 9'2 and J = AlAgg’.
It follows from (6.44) with (2.13), (2.5) and (6.43) that the hyperelastic stresses

are

= %{[l +< (a - 1)] (r2 a?) we} + m
+ 277(1’3‘2'1) (14) (A? cos2 9 + f'2 sin2 Q + Alf, Sin 29) ,

03;,” = gA? [1+ g (A? — 1 + f'2 + g’2)] + h(J),

' It 27 .
03;,” = 7 [1 + g (A? + A? — 1)] 9’2 + h(J) + 7(1)}ib(14)g’2 Sin? 0, (6.52)
0g”: = 0%” = ~3— [1+ 6 (Ag - 1)] f’g’ + 3741), ib(14)g'sinQ (A1 ms!) + f, sin 9) ,
012,53 = 0373...”: = 0,

where h(J) = dH(J)/dJ and (1)/jib (I4) = d<I>fib(I4)/d14 = 14 — 1.

It is to be remarked here that the assumption upon the fluid diffusion direction
in (6.50) can now be verified by the use of (6.16) with &f = 0 by making use of the
fact that p1 and \Ilcm- are independent of planar coordinates a; and y. Namely, since
p1 and \Ilcur in (6.16) vary only in the z-direction, it may then be concluded that the
vector b must be in the z-direction. Therefore it must then follow from (6.31) that
w is similarly aligned in the z-direction.

In view of (6.52) the governing equation (6.30) is automatically satisfied in the

137

y direction, whereas the governing equations in the a: and 2 directions give

d .
Egg” 2 0, (6.53)
and
d - avf
Eagzm + W = 0. (6.54)

It is remarked here that the derivatives d/dz in the Eulerian system are eventu-
ally transformed into derivatives in Lagrangian system d/dZ through the connection
d/dZ = d/dz - dz/dZ = g' - d/dz.

The fluid velocity is now determined by the mass conservation equation (6.4).

For steady-state flow, the integral form of (6.4) in the z direction immediately gives

f
f_i=£flL
v _pf 1—1/J’ (6.55)

where q is an as yet unknown constant which can in turn be related to the fluid mass
flux.
We now turn to consider the boundary conditions. At the bottom (Z = O), the

solid deformation is constrained, leading to
f®=m mm

g®=0 (Mfl

Based on (6.29), the boundary conditions for this steady-state flow problem are
0323(0) = To, 033117? ) = 0, (6-58)

where T0 is the traction reaction due to the support provided by the rigid porous

138

 

plate. It follows from the equilibrium equation (6.12) and the boundary condition
(6.22) that
To = PH — P0. (6.59)

At the top there is no support structure so that t = 0, hence requiring in addition

to (6.58)2 that 0%i$(H) = o and 0%”(H) = o. The first of these, aggxw) = 0, is
automatic by (6.52). It is also to be remarked that (6.59) shows that the boundary
conditions (6.58) are only dependent on the pressure difference PH — P0. Thus the
fluid flow is only affected by the pressure difference PH — P0 rather than the specific
values of fluid pressure PH on the upstream side or PO on the downstream side.

In order to study the problem in more detail we shall temporarily imagine a
broader class of boundary conditions at Z = H than that given by the above condition

0%i3(H) = 0. Narnely we consider either specifying the shear displacement fit with

associated boundary condition

f (H ) = fh (5-60)

or else we may specify the shear traction 1'3 with associated boundary condition
agixw) = ts, (6.61)

In particular, (6.61) with 7'3 = 0 recovers the boundary condition of physical interest
for a free surface

ag'gmw) = o. (6.62)

In principle, (6.53) with (6.52) provides an ordinary differential equation involv-
ing f’(Z), g'(Z), f”(Z) and g”(Z). Substituting from (6.55) into (6.54) and again
using (6.52) provides another such second-order ordinary differential equation for
f (Z) and 9(2). This second ordinary differential equation also contains the as yet

unknown constant q from (6.55). Thus for a well-posed problem we anticipate the

139

need for five boundary conditions. There are indeed 5 boundary conditions, namely
four from (6.56)-(6.58) and one from either (6.60) or (6.61). Therefore this seems to
be a well-posed boundary value problem at least from a mathematical point of view.
Hence it should be possible to set up a numerical routine to solve this problem.

An interesting aspect to this problem is that in the absence of fiber reinforcing
('7 = 0) the solution will involve f = 0. This can be seen by solving (6.53) with
boundary condition (6.56). This f = 0 solution also corresponds to the type of

situation considered by previous researchers in [46] and [71].

6.4.1 Numerical Analysis

We take the following nondimensional forms

__Z -_f __9 __7 —_M
Z-Htf—Htg—H,7—#,AL—H, mm)
and _
f/__d£_fl__l I=i€=£=g1
_dZ dZ ’ dZ dZ ’ (6.64)
”=£1f_’=i=ifn ”=19:_ d9, _1-”

 

fi—Hfl—fiw

With the use of (6.52)3, (6.52)4 and (6.55), the governing equations (6.53) and (6.54)
become two nondimensional second order ordinary differential equations, which are
too long to be displayed here and so given in Appendix A.

A 4-th order Runge-Kutta Scheme and a shooting method are adopted to solve
these two nondimensional governing equations. In particular it is most eflicient to
assign the constant q and to calculate PH -- P0. Some of the material parameters
are obtained from Paul and Ebra—Lima [24]: [.l. = 2.375 x 106 dyne/cmz, X = 0.425,
M = 100p, p5 = 0.869 g/cm3. The initial thickness of the slab is taken to be the
same as that of the experiments described in [24], namely H = 0.0265 cm.

First we consider the case where there are no fibers (7 = 0) and hence no shear

140

displacement (f (Z) ;—= 0). This case corresponds to the experiments conducted in [24].
Since f is formally eliminated in this case, the three relevant boundary conditions are

then given by
9(0) = 0, a$i$<m = PH — Po, oé’é’WH) = 0. (6.65)

Numerical solutions associated with the above boundary conditions are shown in
Figure 6.6a and compared with the experimental data in [24]. Thus the shooting
involves matching to the value 027;” (H) = 0. In this case, since f’ = 0 and 'y = 0 in
(6.52)4, it follows that 0727;,“ = 0 everywhere for this solution.

The values 6 = 0, n = 5, a = 1.556 x 1019 gm/cc - day give a good fit to the
experimental data from [24] (solid curve in Figure 6.6(a)). For the sake of comparison
these values of n and a are then used to obtain the the flux-pressure curve associated
with 6 = 1 (dashed line in Figure 6.6(a)). This comparison shows that the 5 = 1
model predicts far less mass flux at any fixed pressure difference than does the é = 0
neo-Hookean model for the same n and a. The variation of the nondimensional normal
displacement (z — 2) / H along the plate thickness is also shown in Figure 6.7(a) for
E = 0 at the pressure difference PH — P0 = 200psi.

Conversely for .5 = 1 we find that the values of n = 16, a = 3.014 x 1020 gym/cc.
day give a good fit to the same set of experimental data (solid curve in Figure 6.6(b)).
These values of n and a are then used to obtain the flux-pressure curve associated
with 5 = 0 (dashed curve in Figure 6.6(b)). This shows that the E = 0 model predicts
much greater mass flux at any fixed pressure difference than the model of 5 = 1.
The nondimensional normal displacement (z - 2) / H along the thickness direction is
depicted in Figure 6.7(b) for the case where g = 1 and PH — P0 = 200psi.

Figure 6.7 for both 5 = 0 and E = 1 shows that the isotropic layer is contracted in

the thickness direction from its free swelling value. The maximum normal displace-

141

 

_s
N

 

 

 

 

 

 

 

310‘
d
"d
N 8-
E
o
\
3i 6
xo
&
U 4' I .4
g —Numerical Results (5 = 0)
E: 2[ "-Numerical Results (g = 1) q
° Experimental Data
0 200 400

Pressure Difference PH — Po (psi)
(a) 7 = 0, n = 5, a = 1.556 x 1019 gm/cc - day

 

 

 

 

 

 

 

 

 

1 2 a . .
---Numerical Results (5 = 0)
—Numerical Results = 1

310. : E . (é ) a
:3 . o xperimental Data
"‘3 :
mg 8r :
\ 5
c8, 6' :
“so '
g l
a 4- :
El 2_ E
I
0 fl . .
O 200 400

Pressure Difference PH — P0 (psi)
(b) ’y = 0, n =16,a = 3.014 ><1020 gm/cc . day

Figure 6.6: Comparison of the numerical results with the experimental data [24]

142

 

 

 

—PH —: P0=2OO psi

 

 

 

 

 

-0.20

05
Z/H

1

(a) 7 = 0, g = 0: n: 5, a = 1.556 x 1019 gm/cc-day

 

 

 

—PH JP.) = 200psz‘

 

 

 

 

 

O
-0.005-
m -0.01-
>\
“I" -0.015'
N
V -0.02-
-0.025~
-0.03

0

0.5
Z/H

1

(b) '7 = 0, 5 = 1: n =16, a = 3.014 x1020 gm/cc-day

Figure 6.7: The variation of nondimensional displacement in the thickness direction
after free swelling at PH — P0 = 200psz' with 'y = 0 for 5 == 0 and 5 = 1.

143

ment associated with the fluid diffusion (which occurs at the upstream boundary
Z = H) with various pressure difference PH — P0 is shown in Figure 6.8 for both
5 = 0 (with n = 5, a = 1.556 x 1019 gm/cc - day) and for 5 = 1 (with n = 16,
a = 3.014 x 1020 gm/cc - day).

 

 

 

 

’0'250 260 460
PH — P0

Figure 6.8: The normal displacement of the free surface vs. pressure difference for
the material parameters in Figure 6.7(a) with 5 = 0 and Figure 6.7(b) with 5 = 1

6.4.2 Boundary Value Problems

Now we begin to study the effect of fiber reinforcement. This in turn requires
the consideration of the shear displacement f (Z) in the :r-direction. The constitutive
model requires specification of f2 and 7. In what follows we take Q = 1r / 4 and consider
the effects of various 7. In the numerical analysis that follows, the material parameter
5 is taken to be 0. Accordingly we will choose the corresponding values of n = 5,

a = 1.556 x 1019 gm/cc . day.

144

Lateral constraint boundary condition at Z = H

First we discuss the cases in which there is a fiber effect (7 ¢ 0) but with the

surface Z = H held fixed in the lateral direction after free swelling. Define

A2 ’1/2 A2 *1/2
P(Z) = A1 sin!) (1 + —?- tan2 a) — A2 cos 12 (1+ —% cot2 {2) Z
*1 )‘2
(6.66)
so that it follows from (6.41) that P(Z) is the pure free swelling shear displacement
that is associated with the fiber reinforcement. The boundary conditions that we first

wish to consider are given by

M) = 0, 9(0) = 0, 0313(0) = PH — Po.
. (6.67)

f(H) — P(H) = 0, 02;.”(H) = 0.
The fourth order Runge-Kutta scheme with the shooting method is again employed
to solve this boundary value problem. The flux-pressure difference relation obtained
from this procedure is depicted in Figure 6.9(a). It shows that if the fiber reinforcing
parameter 7 is increased then a larger pressure difference is required to sustain the
same amount of fluid flux. The shear stress 0%” at the top surface Z = H that is

necessary to maintain this deformation is then plotted in Figure 6.9(b).

We now turn to consider the fiber deformation, and begin by recalling that the
fibers remains straight lines after free swelling since free swelling is a homogeneous
deformation. However, with the pressure difference and fluid flow, the deformation is
no longer homogeneous and so the fibers now become curved (see Figure 6.10). Their

deformed geometry and slope are shown in Figure 6.11.

145

 

 

—X_X
Q N

 

3
U
'U
N 8'
E
U
E
s, 6‘
"HO
&
c: 4

.....
....
—

 

+7 = 0.0

.3325?" —7 = 0.2
$331" - - ~77 = 0.5
(5252’ ------- 7 = 1.0
..... 7 — 2.0

 

 

 

 

 

0

200

400

Pressure Difference PH —— P0 (psi)
(a) Fluid flux q/pg versus pressure difference PH — P0

 

 

 

 

 

 

 

 

 

0.02
C 5
:1—0.02-
> '\,"'-,“. ~‘.~~
§-0.04- _ ‘ ~~~~~~~~~~
s: *3: 33
b _ _ —7: , I‘.‘ .............
0.06 ___7 = (,5 ,,,,,,,,,,,,
------- 7 = 1.0
_0.08 ,,,,, 7 = 2.0 """~,\
—0.1 . + f
0 200 400

Pressure Difference PH — Po (psi)
(b) 0%”(H) versus pressure difference PH — P0

Figure 6.9: Variation of the fluid flux q/pg and the required reactive shear stress

mix
ozx

(H) with pressure difference PH — P0 for various fiber reinforcement coefficients

7 with 5 = 0, n = 5, a = 1.556 x 1019 gm/cc- day when the boundary conditions are

given by (6.67).

146

 

 

 

V////// ///:>////////////2 M

z . z Seepage
y :> y Direction
A Pressure
Swelling x Driven Seepage K

Figure 6.10: The fibers become curved in the inhomogeneous deformation associated
with pressure driven fluid seepage as depicted in panel (c). This is a schematic
representation.

No shear traction boundary condition at Z = H

In the second case, we consider a free surface boundary condition at the t0p of
the layer Z = H. Then the corresponding boundary conditions are
f (0) = 0, 9(0) = 0, 02’2“”(0) = PH - P0,

. (6.68)
077115111): 0, 0229mm = o

The relation between the flux and the pressure difference as determined from numer-
ical solution of this boundary value problem is depicted in Figure 6.12(a). Again, the
numerical results show that stiffer fibers tend to decrease the fluid flux for the same
pressure difference. The additional lateral displacement at the top surface Z = H
after free swelling f (H) — P(H) is also plotted in Figure 6.12(b).

This then motivates us to study the effect of the initial fiber orientation on the
shear displacement of the free surface. To consider this effect, we first keep the fiber
reinforcing coefficient at 7y = 2.0 and vary the pressure differences PH —- P0 (Figure
6.13(a)). The results show that the maximum shearing displacement of the free
surface occurs at fiber orientation angle 0 > 1r/4. If PH — P0 is small then the angle
9 associated with this maximum is very near 7r / 4. As PH — P0 increases the angle

0 that locates the maximum also increases. In the other way, we could instead keep

147

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0 . 1 i
O 0.2 0.4
:13 f / H
(a) Fiber geometry 2 / H versus a: f / H
5 i
4 h "I."
a: , , .—
4é3“‘ :
33- ’,—":.:“‘:j?
N ........... I"
"C3 ’ ......... "f
2 t ‘‘‘‘‘‘ ".I"‘ .1
1 . .
0 0.2 0.4
13 f / H

(b) The slope of the fibers dz / d2: f versus a: f / H

Figure 6.11: The fiber geometry at various seepage rates with 7 = 2.0, 5 = 0, n = 5
and a = 1.556 x 1019 gm/cc - day when the boundary conditions are given by (6.67).
(23f denotes the projected distance onto x-axis between an arbitrary point on the

fiber and the tip (2:0) of the same fiber.)

148

 

 

 

N # CD
t v .
\
s
s
\
\
\

12
*7 z (,0
3 —7 2 (,2
s 10 ___7 2 0.5
N. ....... 7 : 1.0
g 8 ..... 7 Z 2.0 ‘ ,,,,
\\\ "'ZIN.A:T ____
0 T” . ............
3 """
To "Ii-,3...
Q
\
D"
><
3
Lu

 

 

 

O . .

0 200 400
Pressure Difference PH — Po (psi)

(a) Fluid flux q/pg versus pressure difference PH — P0

 

 

 

 

 

 

0.15 .
*7 = 0.0
—7 = 0.2
_ "-7 = 0.5 ‘‘‘‘
i 0.1 ....... 5: ,0 ,,,,, . ________
A -----— = 2.0 ................
E 7 {,3 ...........
L]. 005 ’Zfiiiili;""
E if";
f: C
-0.05

 

 

 

0 200 400
Pressure Difference PH — Pn (psi)
(b) ( f (H ) — F(H))/H versus pressure difference PH — P0

Figure 6.12: Variation of the fluid flux q/pg and shear displacement f(H) — P(H)
with pressure difference PH — P0 for various fiber reinforcement coefficients 7 with
5 = 0, n = 5, a = 1.556 x 1019 gm/cc - day when the boundary conditions are given
by (6.68).

149

PH — P0 fixed (at 200psi in Figure 6.13(b)) and vary the fiber reinforcing coefficient
7. It is again found that the maximum shear occurs at an angle Q greater than 1r/4.
However it is now found that this angle is close to 77/4 when 7 is large. The angle

increases further beyond 7r /4 as 7 decreases.

Define

2 " 1/2 2 — 1 / 2

T(Z) = A1(1+%cot29) sinQ+A2 (1+§%tan2 Q) 0039 Z,
2 1 (6.69)
where T(Z) is the normal displacement caused by the free swelling deformation (6.41).
Then the dependence of the normal displacement on original fiber orientation for
different pressure difference is shown in Figure 6.14(a) and for different 7 is shown in
Figure 6.14(b).

Again in this case, the fibers become curved and their deformed geometry is
shown in Figure 6.15. In particular, comparing Figure 6.11(a) to Figure 6.l5(a)
shows that the traction free boundary condition gives rise to more overall shortening

of the layer in the z-direction for the same amount of fluid seepage. A corresponding

conclusion holds if the pressure difference is kept fixed.

6.5 General Formulation for Non Steady State Dif-
fusion

In this section, we discuss the time dependent fluid seepage problem in which
case the fluid velocity changes with time. The various field variables are now functions
of (z, t). In particular we consider the following initial-boundary value problem. We
start with a steady state flow problem as described in Section 6.4 with an initial
upstream pressure PH’ an initial downstream pressure P6, and hence an initial fluid

pressure difference PI — PI . We then suddenly change the upstream and downstream
H 0

150

 

0.12~

 

 

 

 

 

 

 

     
     

 

 

 

 

:7 0
PH — P0 -.—_- 0
0 S2 7r/2
(a) 7 = 2.0 with different PH —— P0 = O
0.08
5 7 = 4.0
t: 0.06’ *
>
5 0.04-
Li
I
5
R
V O
-0.02

0 Q 7r/2

(b) PH — P0 = 200psz' with various 7

Figure 6.13: The shear displacements as a function of the initial fiber orientation
under a free surface traction boundary condition at Z = H.

151

 

 

 

 

 

 

 

' 0 fl 7r/2

(a) 7 = 2.0 with different PH — P0

 

-0.02-
m -0.04-
3-006-
5-008

i“ -O.1-
,L-0.12_

 

 

 

 

 

 

 

0 f2 7r/2

(b) PH — P0 = 200psz' with various 7

Figure 6.14: The normal displacements as a function of initial fiber orientations with
a free surface traction boundary condition at Z = H.

152

 

1.5

 

z/If

Ag

(15-

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

315
3 _ ’ ’ l ’ ,,,”, ........................... ;
2 _ 5 ; v ' ' ’ ’ ’ "i”
§ 2. ,,,". _q/pg =0. 0
T3 1 5: "'q/ 105 =30
' " ------- f=60
1’ Q/Po .
1i". """ q/p6=9.0
01%) (f5 1
(Bf / H

(b) The slope of the fibers dz/dxf versus :1: f / H

Figure 6.15: The fiber geometry at various seepage rates with 7 = 2.0, 5 = 0, n = 5
and a = 1.556 x 1019 gm/ cc - day when the boundary conditions are given by (6.68).
The common slope value at Z = H in (b) is a direct consequence of the boundary
conditions (6.68)4 and (6.68)5.

153

pressures to PH and P67 so that the pressure difference abruptly becomes P15 — P5.
The evolution of the solid deformation as well as the fluid flux is of interest. As
t —> 00 one would anticipate that the solution to this problem would approach the
steady state diffusion associated with PH — P6".

The boundary condition for t > 0 follow again from (6.58) and (6.59) as
mist: __ F F mix _
0’22 (0,t)—PH—PO , Uzz (H,t)—0. (6.70)

Note that it makes no difference whether the pressure change AP = (P5 — Pg.) —
(PH — Pg) occurs all on the upstream side (Z = H), all on the downstream side
(Z = 0), or in some combination.

Now the solid velocity v3 is no longer zero but v3 = vgéx + vgéz instead.
Similarly the fluid velocity vf = vgcféx + vzéz We now recall that the dissipative
interactive force 5 as well as seepage velocity w only has the component in the z
direction (see the discussion after (6.52)). Since w = vf —v3 = (v; —f)éx+(v£ —g)éz

it follows that v; — f = 0. Here we denote

6 Z,t 6 Z,t

The mass balance equation (6.4) no longer gives the steady state condition (6.55) but

5% (1 _ f) + 5‘: [(1 — f) 61:] = 0. (6.72)

The stress equilibrium equations now become

instead now gives

52657;.” = 0, (6.73)

154

and

a mix Muir-.51) _
82:0“ + (J ‘1)" — 0. (6.74)

For the convenience of the numerical scheme it is necessary to transform the Eu-
lerian partial derivatives into Lagrangian partial derivatives in the following analysis.
For an arbitrary variable cp = <p(z, t) = 1p[z(Z, t), t], its Lagrangian partial derivatives

could be determined from their Eulerian forms with the use of chain rule:

 

 

 

 

 

 

 

 

 

 

 

 

 

ggzziv 2.53:2 :52; +323 912,0. (g(Z,t)=%a@[ ).
Therefore it follows that
%.=S_§JZ%T> (6.75)
and

 

 

where the second equality of (6.76) has made use of (6.75).

6.5. 1 Isotropic Layer

In this section, we consider the transient diffusion problem described in Section
6.5 for the special case in which the layer is isotropic (7 = 0). Thus A1 = A2 =
A = A1 and there is again no shear displacement, i.e., f (Z) 5 0. Thus the stress
equilibrium equation (6.73) in the 117-direction is trivial. Making use of (6.75) to

evaluate 8022 235/82: it is found that the governing equation (6. 74) in the z- -direction

155

 

gives

 

A4 A2 1 n
1+ {(212 — 1) + —h’(A29’> g” + a ( ) (v! — 1119’: o, (6.77)
u u A
where h’(J) = dh/dJ.

Using (6.75) and (6.76) to express (6.72) in terms of Lagrangian partial deriva-
tives, and then eliminating 72'; between that result and (6.72) gives a partial differ-
ential equation for g(Z , t) that is third order in space and first order in time. This
equation is of the form
1

CF2(9’)9”2, (6.78)

.I 1 I III
=—F
g c 1(g)g +

where F1(g') and F2 (9’) are both unit free functions of the single field variable g’, and
C = 01A4 / [1. Note that C has units of sec/cm2. These functions are given in Appendix
B The most obvious boundary condition is given by the fixed displacement boundary
condition g(O, t) E 0 at the downstream location Z = 0. It is also possible to obtain
boundary values for g’(0, t) and g’ (H, t). These are found by solving equation (6.70)
using (6.52) with 7 = 0 and J = Azg’. This is done numerically and does not present
any special difficulty. It is also necessary to assign an initial condition for g(Z, 0) for
0 < Z < H. This is given by the steady state solution associated with the pressure
difference P}! — POI. With respect to the (z, t) domain there is now a discontinuity
in 9’ (Z, t) at the corner (Z, t) = (0,0) because t11_nq)g’(0, t) comes from the boundary
condition while 2111110 g’(Z, 0) comes from the initial condition. Note also that no such
discontinuity in 9' (Z, t) occurs at the corner (Z, t) = (H, 0).

We partition the space domain by using a mesh Z0,...,Z I and time domain using
a mesh t0,....,tN. In particular, at each time step there are I — 1 internal nodes
(i = 1, 2, ..., I — 1) since i = 0 and i = I are boundary nodes. We assume a uniform

partition both in space and in time, so the difference between two consecutive space

156

 

 

Boundary

 

condition

\ of specified

8 UN)

 

 

 

 

 

 

 

 

 

n+1

n a“ .‘ l.‘ 1"“.

 

 

 

 
    

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Boundary /
condition
of specified ‘_ Initial.
g(0,t) and 2 condition
g ‘(0 0 matches

’ 1 boundary

\r 47/ condition
0

1 2 i-l i i+l 1+2 I-lI

Initial Condition of specified g(Z ,0)

 

 

 

 

Figure 6.16: An explicit forward time-stepping scheme

points will be AZ and between two consecutive time points will be At. The points
g(Z,-,tn) =9? (i=0,1,2,...,I, n=0,1,2,...N)

will represent the numerical approximation of g(Z, t). Using the initial and boundary
condition the numerical solution may proceed by an explicit forward time-stepping
scheme using a 6 node stencil to discretize (6.78) as shown in Figure 6.16. In particular
suppose the stencil involves two nodes (i, n + 1), (i + 1, n + 1) on the time step 71. + 1
and four nodes (i — 1,n), (i,n), (i + 1,71), (2' + 2,71.) on the time step 72.. These six
nodes comprises a stencil, whose ”center” is at (i + %, n + %) as represented by the
stars in Figure 6.16. The discrete form of (6.78) on the basis of this stencil is now

written as
:left =
i+%,n+%

[I]

”'9?“ 1 (6.79)
i+2,n+2

with

157

 

 

 

1 n+1 1 .+1
:left = H(gz'fi ‘ 9&1) ‘ Em? ‘ 9?)
i+%,n+% AZ
:Tight 2 F1 9ft+1 ‘ 9? 92n+2 " 392n+1 + 39f - gen—1 (6 80)
245.1%? AZ C(AZ)3 ‘

1 9"“ - g" 1 2
Z Z
+ 'C-F2 (T) [W(gan _ 910+1 — 9? + 9:11) ~

Assume that numerical values have been obtained for g(Z, t) at time step 71 and
observe that gg'H is known from the boundary condition. It is then desired to obtain
numerical values for g(Z, t) at time step 71 + 1 by using the explicit scheme (6.79) for
i=0, 1, 2, ..., I-1. It is to be noted that the use of (6.79) for both i = 0 and i = I — 1
requires the consideration of a virtual node as shown by the dashed triangles in Figure

6.16. The values 9’11 and g?+1 of these two virtual nodes are determined from the
known boundary values g’(0, t) and g’ (H, t) through

I 1 n n I _ 1 n
g (0.t) = $7191 -9_1). and 9 (Hi) — 55-2-6?“ -91_1)- (6-81)

Numerical solution of the partial differential equation (6.78) is then obtained
through the forward time stepping scheme (6.79). In particular, at the end of time
step 71. the values of g? are known for i = 0, 1, 2, ..., I. The values 931 and g?+1 then
follow from (6.81). The value gg+1 is also known from the boundary condition. Then
(6.79) fori = 0,1,2, ..., I—1 gives I equations for the I unknowns 9271—51 (i = 1, 2, ..., I).
These equations are easily solved in a sequential manner since the only unknown in
(6.79) for i = 0 is 9’1“”1, the only new unknown in (6.79) for i = 1 is then 93+1, etc.

Figure 6.17 shows the results of such a numerical calculation. The initial pressure

difference in this example is PH — P6 = 50psi. The pressure difference is then

suddenly changed to P}; — POP = 200psi. The approach to the steady state solution

158

 

0 .

 

 

 

 

\‘~.\Initia/l Steady State
‘~\ 2 _ —8
-0.04- ‘~~~t./£§'{{)— “0 .
m ......... .
>
E -0.08 2 X 10—7 ‘
I
3 . 5 x 10'7
_0'12' 1 x 10-6
Final Steady State _ >5 10"
'0'160 02 0:4 0:6 08 1

Z/H

Figure 6.17: Change in compaction during transient diffusion in the absence of fibers
(7 = 0) for a material with 5 = 0, M / a = 100. In this simulation the steady state
seepage with PH -— P0 = 50psi is abruptly altered to PH — P0 = 200psi.

159

is evident.
Base on the results in Figure 6.17, the profile of fluid density ratio pf could be
determined as in Figure 6.18. In this case pf = (1 — J_1) pf = (1— A—2g’_1)p5.

 

Initial Steady State PI}; — Pop 2 50 psi

0.8-

   
 
  

Final Steady State
P}; — Pop 2 200 psi _

t/(CHz) = 5 x 10‘8,2 x 10‘7,
5 x 10-7, 1 x 10-6, 5 x 10-6.

 

.1

 

0 0.2 014 06 0.8 1
Z/H

Figure 6.18: The variation of pf/pg along the thickness for the simulation considered
in Figure 6.17.

The variation of normal stresses egg”: and egg,“ along the thickness direction is

shown in (Figure 6.19). Note here the symmetry gives 0%” = 05,7393.

The total mass of fluid per unit current cross-section area within the layer 2 can
also be obtained by 2 = fg‘ pfdz = [3" pfg’dZ = fol pfg'HdZ (Figure 6.20). The
numerical routine shows good convergence of the function g to its large time steady
state value.

The example connected to Figure 6.17 — Figure 6.20 concerns the transition from
one steady state flow to another steady state flow. Alternatively, one can begin with

an initial state that involves free swelling and hence no flow by virtue of PH — P0 = 0.

Then by abruptly raising the value of PH — P0 one obtains a transient flow condition

160

 

   
  

 

 

 

 

   
   

 

Initial Steady State Pg — P6? = 50 psi
0* /
i _____ ..
\ "0.1 Final Steady State
g 8 Pg," — P5“ = 200 psi
g—oz
_0 3 t/(C’H2) = 5 x 10‘s, 2 x 10‘7,
' 5 x 10-7,1 x 10-6,5 x 10-6.
—0.4 . . . .
0 0.2 0.4 0.6 0.8 1
Z/H
(a) agix/a versus Z/H
0 Initial Steady State P}; — Pf = 50 psi
3. """"
8\ “0.2 ‘ Final Steady State‘
E § P}; — Pf, = 200 psi
b
_0'4 t/(CH2) = 5 x 10‘s, 2 x 10”, l
5 x 104,1 x 10-6,5 x 10-6.
-0.6

 

 

0 0:2 04 0.6 018 1
Z/H
(b) agnzfx/a versus Z/H

Figure 6.19: The variation of normal stresses along the thickness direction for the
simulation considered in Figure 6.17.

161

 

1.36
1.34 \

Initial Steady State PI]; — Pf = 50 psi
A 1.32 .

m 1.3

Kg? ‘
< 1,28-A Final Steady State
W P}; — P5 = 200 psi

1.26b

1.24

1“220 0:5 1 1:5 2 25 3
t/(CHz) x10-5

Figure 6.20: The amount of fluid in the layer changes with time: from its initial
steady state value (where PH — P0 = 50psi), and finally approaching the final steady
state value (where PH — P0 = 200psi).

 

...........................

 

 

 

162

that will approach a steady state for large time. Figure 6.21 shows this phenomena

in a manner similar to Figure 6.20 for a variety of P11; — POF values.

1.4

 

    
 
 
  
  
 

Steady {State Value for free I
swelling (P11; — Pf = 0)
1.35 ........ Steady State Value for P5 — Po“ = 50 psi

 

 

A Steady State Value for P5 — Pf = 100 psi
m 1 3 ......... .7
“so
Q.
V
a 1 25 Steady State Value for P5 — P6? = 200 psi _

 

Steady State Value for P}; - POF = 300 psi

 

 

 

 

 

1 2 ,_ ..............
Steady State Value for P}; — P0F = 400 psi
1'1550 0:4 0:8 1:2 1:6 2

t/(CHZ) x10-5

Figure 6.21: The change of interstitial fluid amount with time for various pressure
difference. Free swelling state (PH — P0 = 0) is abruptly altered when PH — P0 is
suddenly changed to 50psi, 100psi, 200psi, 300psi, and 400psi, respectively.

6.5.2 Fiber Reinforced Layer

Now we consider the transient diffusion problem for a fiber reinforced layer. As
in the previous section, we suddenly increase the pressure difference from P}! — P6
to P5 — P6? and analyze the evolution of the solid deformation and the fluid flux.

The stress equilibrium equations are the same as in (6.73) and (6.74). The fluid
mass balance equation is again given by (6.72). By substituting from (6.74) into
(6.72), the evolution equation regarding the time derivative of g(Z, t) is obtained as
follows:

, 1
g’ = 51730“, f”, f”’,g’,g”,g”’), (6.82)

163

where C' = aA¥/\%/p, and F3 is a complicated function of f’, f”, f’”,g’,g” and 9’”
that is too cumbersome to display here. It is given in Appendix C. This equation is
the generalization of the previous equation (6.78) to the case of fiber reinforcement.

It is immediate that one set of boundary conditions is
f(0,t) = 0, and g(O, t) = 0. (6.83)

The remaining traction boundary conditions are as follows. First condition (6.70)

continues to hold, namely

 

agixm, t) = pg — POF, 633w, t) = o. (6.84) 1

The remaining traction boundary conditions follow from the requirement that a“ (H, t) =

0 which, in conjunction with (6.53), gives
agi$(0,t) = 0, agi$(H,t) = 0. (6.85)
In fact (6.53) gives the immediate first integral
”g(Z, t) = 0. (6.86)

Note also that by using (6.52) it follows that (6.84)1 and (6.85)1 provide two equa-
tions for f’ (0, t) and g’(0, t). Similarly (6.84)2 and (6.85)2 provide two equations for
f’ (H, t) and g’ (H, t). These equations can be solved numerically for constant values
f’(0,t), g’(0, t), f’ (H, t) and g’ (H, t) prior to the consideration of the time-stepping
numerical routine.

For 0 < Z < H, initial values of f (Z, t) and g(Z, t) are given by the steady state

solution to the problem with fiber reinforcement. With these boundary and initial

164

conditions, we may proceed with the numerical iteration by using the same forward
time stepping scheme as in Figure 6.16. Observe that the only time derivatives that
appear in the formulation are 9 and g’, and these time derivatives only appear in
equation (6.82). Therefore, we obtain nodal values gy‘fl for i = 1, 2, ..., I — 1, I when
we move from time step ii to time step n + 1. It is then necessary to use (6.86) to
solve for fin+1 for i = 1,2, ...,I — 1,] from these gzm—l. Similar virtual nodes for
i = —1 and i = I + 1 are then created in the time-stepping stencil scheme. Namely,
the values of fill”, 9711-1, fill—11 and 95111 will be introduced to complete the scheme
with the use of the boundary values of f’(0,t), f’(H, t), g'(0, t) and g'(H, t) as in
(6.81).

Numerical approximation of the profiles of f (Z, t) and g(Z, t) is shown in Figure
6.22. The parameters are chosen as follows: 5 = 0, the fiber reinforcement coeflicient
7) = 1.0 and the original orientation angle of the fibers 9 = 1r/4. Recalling the
definitions in (6.66) and (6.69), the function P(Z) and T(Z) represent, respectively,
the pure shear and normal deformation caused by the fiber-reinforced free swelling.
Figure 6.22(a) and Figure 6.22(b) then depict the evolution of f (Z) — P(Z) (shear
displacement due to pure pressure difference) and 9(2) — T(Z) (contraction due to
pure pressure difference), respectively.

Comparison of Figure 6.22 with Figure 6.17 indicates that the presence of fibers
causes the time for the system to achieve its final steady state to increase by nearly
two orders of magnitude over that of the isotropic case. The shear displacement and
the normal contraction for the fiber-reinforced transient diffusion process with an
increased fiber coefficient 7) = 10 is shown in Figure 6.23. The overall way in which
the transient deformation approaches the steady state deformation is found to be
relatively insensitive to the the initial fiber orientation Q as seen in Figures 6.24 and

6.25.

165

0.06

 

 

 

 

 

 

 

 

 

 

 

 

 

0 05 Figal Stgady State 5 x 10:: - .. -.
. " P —P =200psi ,x” 4
I H 0 \ x” 2 x 10‘4
$0.04. I,” 1x10—4‘
:003_ ” 5X10-5
,L ' t/(C’Hz) = 1 x10"5
5, 0.02~
b , —————————— 7"'

0.01 ' '''''' Initial Steady State‘

0 ,x'f . Pg—ijztilOpsi
O 0.2 0.4 0.6 0.8 1
Z/H
(a) (f(Z) — F(Z))/H versus Z/H
0 “~\ ' ' Initiai Steadgj State
Pg —Pg“ = 50 psi

I -o.02- .........
E: t (0112)";13516-75
N -0.04-
g: . 5 x 10‘5
ii? -0.06' ‘ \ 1 x 10-4.
é ‘~\ 2 x 10-4

—0-08" Final Steady Statefx‘~ 5 x 10—4"

P}; — P67 = 200 psi ‘ ________
-0 1

'0 0:2 0:4 0:6 0:8 1
Z/H

(b) (g(Z) — T(Z))/H versus Z/H

Figure 6.22: The shear and normal contraction associated with fiber reinforced tran-

sient diffusion where '7 = 1.0, Q = 45°, and E = 0.

166

 

0.08

 

 

 

 

 

 

 

Final Steady State 5 X 10:4: ’’’’’’’’ ,
Pf; — P5” = 200 psi 2 x 10-4
0-06' \/ . 1 x 10-4‘
’2’: I 5 x 10‘5
Di
5: 0'04" t/(C’Hz) = 1 x 10-5‘
'2‘,
0.02- ___________________
""""""" I nitial Steady State
P5-P{=50 psi
00 0.2 0.4 H 0.6 0.8 1
(a) (f(Z)— I‘(Z))/H versus Z/H
0 “‘~
-o.o1- _____________
E -0.02-
S -o.03-
T
g; -o.04~
-0.05" ‘~“‘
—o.06- ____________
‘0'070 0:2 0:4 0:6 0:8 1
Z/H

(b) (g(Z) — T(Z))/H versus Z/H

Figure 6.23: The shear and normal contraction associated with fiber reinforced tran-
sient diffusion where ’7 = 10.0, (2 = 45°, and g = 0.

167

 

0.015

 

 

 

 

 

 

 

; 0.01 .
S /
[— I
I I
K; I
g 0.005- ’ _ _‘
0 I i . . .
O 0.2 0.4 0.6 0.8 1
Z/H
(a) (f(Z) — P(Z))/H versus Z/H for '7 =1.0 and S2 = 15°
0 \
-o.02- ‘ r . ‘
5-004» ~ ‘~~-___)
Q \
a -0.06 x
a” \
g —o.08~
9 —o.1 -
—o.12» \ . ~
-014 . . . .
0 0.2 0.4 0.6 0.8 1
Z/H

(b) (g(Z) — T(Z))/H versus Z/H for 7y = 1.0 and Q = 15°

168

 

0.04 .

 

 

 

 

 

 

 

0.03’
I
8
ET 0.02'
S:
V 0.01? _____
o ’ ’ .
0 0.5 1
Z/H
(c) (f(Z) — I‘(Z))/H versus Z/H for '7 =1.0 and Q = 75°
0 s
-o.02- ‘ ‘ ----
I
8:
>7 —o.04r
S
9
-0.06"
-0.08 r
O 0.5 1
Z/H

(d) (g(Z) — T(Z))/H versus Z/H for '7 =1.0 and Q = 75°

Figure 6.24: The shear and normal contraction associated with fiber reinforced tran-
sient diffusion ’where '7 = 1.0, Q = 15° and f2 = 75°

169

Z)-F(Z))/H

v
I.—
v

(a) (f(Z) — F(Z))/H versus Z/H for 7y = 1.0 and Q = 40°

(g(Z)"Y(Z))/H

(b) (g(Z) — T(Z))/H versus Z/H for '7 =1.0 and Q = 40°

 

0.06

 

 

 

0.05- , , , — - -
0.04- x”
003. z”
0.02
OD1~ ,s”””
0o ’ 0:5 1
Z/H

0
-0.02'
-0.04-
-0.06-
-0.08*

-0.1 -

 

 

‘~
~
~~
---

 

 

-0.12
0

170

0:5
Z/H

1

 

 

0.06

 

 

 

 

 

 

 

 

0.05- ,z ’
’l
3.; 0.04- ’
E
L... 0.03-
K? .
g 0.02. » 4
0.01 " xx"
0 ’ . 9
0 0.5 1 ~
ZIH
(c) (f(Z) — F(Z))/H versus Z/H for 7y = 1.0 and Q = 50°
0 \ -
—0.02' ‘ ~~~~~~~ J
I - -
§ -o.04-
*r
g -o.oai
92 x
-0.08- ‘~. ‘
-o.1 1
o 0.5 1
Z/H

(d) (g(Z) — T(Z))/H versus Z/H for '7 21.0 and f2 = 50°

Figure 6.25: The shear and normal contraction associated with fiber reinforced tran-
sient diffusion where '7 = 1.0, Q = 40° and Q = 50°

171

 

6.6 Summary

In this chapter, We study the pressure driven fluid diffusion through a hyperelas-
tic porous layer. Both steady state and transient diffusion cases are considered with

particular interests in the effects of fiber reinforcement on the diffusion process.

0 It is found that the original fiber orientation angle 0 is changed to to during the
homogeneous deformation of free swelling. This angle change w—Q is monotonic
with the fiber reinforcing coefficient '7 for a certain angle 9 (0 < 9 < 7r / 2). For
a fixed fiber reinforcing coeflicient 7, the angle change exhibits a maximum

value (greater than 7r/4 if '7 > 0) for different values of Q.

o For steady state pressure driven fluid seepage in the absence of fiber reinforcing
it is found that the model (6.52) with g = 0 providas a good match to experi-
mental curve in [24] by taking n = 5 and a = 1.556 x lolggm/cc - day. These

values are then used in the fiber reinforced model (7 aé 0).

o For steady state pressure driven fluid seepage in the fiber reinforced material,
the resulting nonhomogeneous deformation causes the fibers to be curved. Both
lateral constraint boundary condition and no shear traction boundary condition
at the surface Z = H are considered. In the first case, it is found as predicted
that a) both flux and the shearing stress at Z = H increases monotonically
with the pressure difference; b) the fiber reinforcement causes more fluid flux
and greater shearing stress at Z = H with a given pressure difference; 0) the
fibers get further curved as 7 increases. In the second case, it is again found that
a) both flux and the shearing displacement at Z = H increases monotonically
with the pressure difference; b) the fiber reinforcement causes more fluid flux
and greater shearing displacement at Z = H with a certain pressure difference;
c) the fibers get further curved as 77 increases but with a common slope at

Z = H. The effects of both pressure difference PH — P0 and 5' on the shearing

172

displacement at Z = H are studied.

0 Transient diffusion process is studied in the problem where a abrupt pressure
difference change is applied. Numerical solutions are obtained by using an
explicit forward time-stepping scheme. It has shown that the presence of fibers
causes the time for the system to achieve its final steady state to increase by

nearly two orders of magnitude over that of the isotropic case.

173

Chapter 7

Conclusion

In this dissertation, we focused on theoretical modeling and numerical simula-
tion on the mechanical response of swollen elastomeric gels. To do this, we first
introduce a hyperelastic constitutive theory involving the generalized Flory-Huggins
framework. We then consider the static equilibrium response including homogeneous
and inhomogeneous deformation induced by different mechanical loading conditions.
In particular we consider the equilibrium states in which these mechanical loadings
can lead to both fluid loss (swelling reduction) and fluid gain (swelling increase). For
loadings that give fluid gain, we consider the transition from the state of liquid sat-
uration to that of nonsaturation at a particular point on the quasi-static load path.
Numerical simulations are introduced to solve the stress equilibrium equations. In

what follows list the major numerical results.

0 In Chapter 3, the homogeneous deformation is first discussed. It is found that
both the stress and the overall volume of the gel are monotone with the stretch
for equibiaxial loading, while the volume change exhibits a local maximum at
a certain uniaxial stretch for the case E = 1. For equitriaxial loading, the
stress is found to be monotonically increasing if g is greater than a critical value

Ecrita and exhibits a stress maximum if 5 < gem. In all cases, nonsaturated

174

bifurcation always leads to a stiffer response.

In Chapter 4, the inhomogeneous deformation of everting a swollen tube subject
to axial loading is studies. Numerical results show that the deformed inner and
outer radii decrease as the axial stretch increases. The overall volume of the
tube is monotone in the axial stretch if é = 0 but exhibits a local maximum as
long as is > 0. The transition to a nonsaturated state preserves the overall gel
volume but alters the distribution of the fluid content within the gel. For the
cases 5 = 0 and 5 = 1 this results in opposite directions of the fluid quasi-static

migration in the everted cylinder.

In Chapter 5, the problem of twisting a swollen tube is considered. Two sorts
of boundary value problems leading to the change in the overall gel volume are
solved. Either the inner or outer surface is fixed radially while let the other
lateral surface be free of normal traction. The case where the inner surface is
free of normal traction gives rise to fluid uptake, which makes the transition to
a nonsaturated state possible. The result show that relation between the twist

angle and the twist moment abruptly stiffers after loss of saturation.

As the second part of this work, the pressure driven fluid diffusion through a

fiber reinforced hyperelastic media is investigated in Chapter 6. To study the relative

motion between the fluid and polymer matrix, we invoke a large deformation bipha-

sic mixture theory that specially deals with separate mechanical balance principles

for each . Both steady-state and non-steady—state diffusion processes are considered.

Numerical analysis involving a forward time—stepping scheme is utilized in order to

obtain time dependent swelling behavior of elastomeric matrix as well as the redistri-

bution of the fluid concentration during diffusion. Numerical simulation shows: a) the

fiber orientation angle is changed during the homogeneous of free swelling; b) steady

state seepage causes the original straight fibers to be curved; c) fiber reinforcement

175

causes more flux and greater deformation of the polymer network; d) the study of
transient diffusion process shows that the presence of fibers causes the convergence
time to increase greatly over that of the isotropic case.

There are several potential directions to expand this work. Firstly, the current
numerical simulations take the values of the material parameters M, x, and ,u directly
from the experimental data in [24] for the specific combination of toluene and vul-
canized rubber. In order to compare with other experimental data on different gels,
adjustment might be made according to the corresponding materials. Furthermore,
this work takes Mooney-Rivlin elastic energy function for the elastic deformation of
the polymer network, and Flory-Huggins’ theory to interpret the mixing entropic and
enthalpic effects between the elastomer and the liquid. However, we are also able to
embed more general material models into the current framework.

Another possible improvement is to consider more complicate geometries, and
also other constitutive models involving viscoelastic behaviors, which many sorts of
polymers bear. To do this, finite element methods might be adopted to solve the

more complicate time and loading rate dependent equations.

176

APPENDIX

177

Appendix A

Nondimensional Form of
Governing Equations of Steady
State Diffusion Problem

By introducing the nondimentionized process in (6.63) and (6.64), the stress

equilibrium equations (6.53) and (6.53) then become

T10“, i’)f” + F202 i’)g7’ = 0 (Al)

and
aqH A%/\§§'2

73(f’,§')97’ + f4(f’.§’)f” + =
Mp5 (4M2? "1)"+1

 

0, (A2)

where f1(f’,§') = 1+§()\%-1)
+ 2’7 [sin4 {2(3f’2 + (7'2) + 3A? sin2 9 sin 2Qf_’ + sin2 9 (3A? cos2 (2 —1)],
1326’, g’) = 25' (A? sin 20 + 25in2 Qf’) sin2 Qg’,
fsa‘as’) = 1 +6 (A? + A; — 1) + Aiigh’mliis’w
+ 2’? [sin4 {2(f’2 + 3572) + A? sin 29 sin2 Qf’ + sin2 Q (A? cos2 {2 — 1)],
$407, 6’) = 27 (A? sin 29 + 2sin2 of) sin2 96', and h’(J) = dh(J)/dJ.

178

 

 

Appendix B

Functions F1(g') and F2(g’)

= [A +Bhl()\29I)](Azgl _1)n+1

F1(g,) 9,2

 

(13.1)

, Bh”(A2g’)A2(,\Zg’ — 1)("+1> + (n + 1)[A + Bh’(A2g’)](A29’ — 1m?
F2(g ) = 9’2

A + Bh’(/\Zg’)]()\2g’ — 1)n+1
9’3 ’

 

_21

 

(B.2)
where A =1+£(2)\2 — 1) and B = All/[1.

179

Appendix C

Function F3(f,, f”, fl”,g,,g”,g,”)

F3(f', f”, f’", 9', g”, 9'”) = WAA + BB) + 00 * DD * 9”

_ n+1
+1—2—
9

(0.1)
(DDIQII +DDgIII),

where

_ n III__ _ n+1 II
AA=(n+1)(J 1) AME/929 (J 1) 9 (23in2nf’+A1sin2o)sin2nf”,

_ n+1
BB = £J_;),_ Sin? 9 [2 31112 9f”2 + (2 sin2 Qf’ + A1 sin 29)f”'],
CC _ (n+1)(J-1)nA1/\2g”g’—2(J-1)n+lgll
_ 9/3 ,

DD = 1 + g(A’f + Ag — 1) + A§A§h’(.z)

 

 

+ 2i‘ysin2 Q[sin2 (2(f’2 + 39'2) + A? cos2 Q — 1 + A1 sin2f2f’],
DD’ = A§A3h"(.1)g” + 47 [ sin4 9( f’ f” + 3g’g”) + A1 sin 29 sin2 of”],
J = AlAgg’ and h”(J) = d2h(J)/dJ2.

180

 

BIBLIOGRAPHY

181

Bibliography

[1] B. Baroli, “Hydrogels for tissue engineering and delievery of tissue-inducing sub-
stances,” J. Pharm. Sci, vol. 96, no. 9, pp. 2197—2223, 2007.

[2] Y. Qiu and K. Park, “Environment-sensitive hydrogels for drug delivery,” Adv.
Drug Deliv. Rev., vol. 53, no. 321-339, 2001.

[3] J. C. Wheeler, J. A. Woods, and M. J. Cox et al., “Evolution of hydrogel polymers
as contact lenses, surface coatings, dressings, and drug delivery systems,” J. Long
Term Efi“. Med. Implants, vol. 6, no. 3—4, pp. 207—217, 1996.

[4] T. Noguchi, T. Yamamuro, M. Oka, P. Kumar, Y. Kotoura, S. H. Hyon,
and Y. Ikada, “Poly(vinyl alcohol) hydrogel as an artificial articular-cartilage-
evaluation of biocompatibility,” J. Appl. Biomater., vol. 2, no. 2, pp. 101-107,
1991.

[5] D. J. Beebe, J. Moore, J. M. Bauer, Q. Yu, R. H. Liu, C. Devadoss, and B. H. Jo,
“Functional hydrogel structures for autonomous flow control inside micro-fluidic
channels,” Nature, vol. 404, 2000.

[6] L. Bromberg and E. Ron, “Temperature-responsive gels and thermogelling poly-
mer matrices for protein and peptide delivery,” Adv. Drug Deliv. Rev, vol. 31,
pp. 197—221, 1998.

[7] D. Schmaljohann, “Thermo— and ph-responsive polymers in drug delivery,” Adv.
Drug Deliv. Rev., vol. 58, pp. 1655—1670, 2006.

[8] E. S. Gil and S. M. Hudson, “Stimuli-responsive polymers and their bioconju-
gates,” Prog. Polym. Sci, vol. 29, pp. 1173—1222, 2004.

[9] G. Gerlach, M. Guenther, J. Sorber, G. Suchaneck, K. F. Arndt, and A. Richter,
“Chemical and ph sensors based on the swelling behavior of hydrogels,” Sens.
Actuators B Chem, vol. 111, pp. 555-561, Nov. 2005.

[10] A. Mamada, T. Tanaka, D. Kungwatchakun, and M. Irie, “Photoinduced phase
transition of gels,” Macromolecules, vol. 23, pp. 1517—1519, 1990.

[11] A. Suzuki, T. Ishii, and Y. Maruyama, “Optical switching in polymer gels,” J.
Appl. Phys, vol. 80, pp. 131—136, 1996.

182

[12] S. K. De, N. R. Aluru, B. Johnson, W. C. Crone, D. J. Beebe, and J. Moore,
“Equilibrium swelling and kinetics of ph-rasponsive hydrogel: models, experi-
ments, and simulations,” J. Microelectromech. Syst, vol. 11, pp. 544—555, Octo—
ber 2002.

[13] L. R. G. Treloar, The physics of rubber elasticity. Oxford Univ. Press, 1975.

[14] G. Gee, “The interaction between rubber and 1iquids.9. the elastic behavior of
dry and swollen rubbers,” Trans. Faraday Soc., vol. 42, no. 8, pp. 585-598, 1946.

[15] S. M. Gumbrell, L. Mullins, and R. S. Rivlin, “Departures of the elastic behavior
of rubbers in simple extension from the kinetic theory,” Trans. Faraday Soc.,
vol. 49, no. 12, pp. 1495-1505, 1953.

[16] L. Mullins, “Determination of degree of crosslinking in natural subber vulcan-
izates: Part iv. stress-strain behavior at large extensions,” J. Appl. Polym. Sci,
vol. 2, no. 6, pp. 257—263, 1959.

[17] G. B. Mckenna, K. M. nynn, and Y. Chen, “Experiments on the elasticity of
dry and swollen networks: implications for the frenkel—flory-rehner hypothesis,”
Macromolecules, vol. 22, no. 4507-4512, 1989.

[18] R. A. Paxton and A. M. Al-Jumaily, “Swelling of polyelectrolyte hydrogels using
a finite element model,” Polymer, vol. 47, pp. 5997—6003, 2006.

[19] O. C. Zienkiewicz, A. H. C. Chan, M. Pastor, D. K. Paul, and T. Shiomi, “Static
and dynamic behaviour of soils: a rational approach to quantitative solutions. i.
fully saturated problems,” Proceedings of the Royal Society of London. Series A,
Mathematical and Physical Sciences, vol. 429, no. 1877, pp. 285—309, 1990.

[20] S. Prager, “Diffusion of hydrocarbons in polyisobutylene,” J. Am. Chem. Soc.,
vol. 73, pp. 4072—4075, 1951.

[21] R. E. Pattle, “Diffusion of 1-benzeneazo—2—naphthol in swollen rubber,” Trans.
Faraday Soc., vol. 62, pp. 1776—1784, 1966.

[22] N. A. Pepper and C. T. Reinhart, “Solute diffusion in swollen membranes 1. a
new theory,” J. Membr. Sci, vol. 15, no. 3, pp. 275—287, 1983.

[23] C. J. Guo, D. D. Kee, and B. Harrison, “Diffusion of organic solvents in rubber
membranes measured via a new permeation cell,” J. Appl. Polym. Sci, vol. 56,
no. 7, pp. 823-829, 2003.

[24] D. R. Paul and O. M. Ebra—Lima, “Pressure-induced diffusion of organic liq-
uids through highly swollen polymer membranes,” J. Appl. Polym. Sci, vol. 14,
no. 2201-2224, 1970.

[25] M. Mooney, “A theory of large elastic deformation,” J. Appl. Phys, vol. 11,
no. 9, pp. 582—592, 1940.

183

[26] R. S. Rivlin, “Large elastic deformations of isotropic materialsl. fundamental
concepts,” Phil. Trans. R. Soc. Lond. A-Math. Phys. Sci, vol. 240, no. 822,
pp. 459—508, 1948.

[27] P. J. Flory, “Theory of elasticity of polymer networks-effect of local constraints
on junctions,” J. Chem. Phys, vol. 66, no. 12, pp. 5720—5729, 1977.

[28] P. J. Flory and B. Erman, “Theory of elasticity of polymer networks. 3.,” Macm-
molecules, vol. 15, no. 3, pp. 800—806, 1982.

[29] R. J. Gaylord and J. F. Douglas, “Rubber elasticity-a scaling approach,” Polym.
Bullet, vol. 18, no. 4, pp. 347—354, 1987.

[30] R. J. Gaylord and J. F. Douglas, “The localization model of rubbber elasticity
.2.,” Polym. Bullet, vol. 23, no. 5, pp. 529—533, 1990.

[31] E. A. Dimarzio, “Contribution to a liquid-like theory of rubber elasticity,” J.
Chem. Phys, vol. 36, no. 6, pp. 1562—1570, 1962.

[32] E. A. Dimarzio, “Contribution to a liquid-like theory of rubber elasticity .2.
existence of a (Ax/\y + AyAz + Azkx) term,” Polymer, vol. 35, no. 9, pp. 1819—
1826, 1994.

[33] E. M. Arruda and M. C. Boyce, “A 3-dimensional constitutive model for the
large stretch behavior of rubber elastic materials,” J. Mech. Phys. Solid, vol. 41,
no. 2, pp. 389—412, 1993.

[34] W. H. Han, F. Horkay, and G. B. Mckenna, “Mechanical and swelling behaviors
of rubber: a comparion of some molecular models with experiment,” Math. M ech.
Solids, vol. 4, pp. 139—167, 1999.

[35] R. S. Rivlin and D. W. Saunders, “Large elastic deformations of isotropic mate-
rials. 7. experiments on the deformation of rubber,” Phil. Trans. R. Soc. Lond.
A-Math. Phys. Sci, vol. 243, no. 865, pp. 251—288, 1951.

[36] H. Pak and P. J. Flory, “Relationship of stress to uniaxial strain in crosslinked
poly(dimethylsiloxane) over the full range from large compression to high elon-
gations,” J. Polym. Sci. B-Polym. Sci, vol. 17, no. 11, pp. 1845—1854, 1979.

[37] G. Allen, M. J. Kirkham, J. Padget, and C. Price, “Thermodynamics of rubber
elasticity at constant volume,” Trans. Faraday Soc., vol. 67 , no. 581, pp. 1278—
1292, 1971.

[38] P. J. Flory and Y. I. Tatara, “Elastic free-energy and elastic equation of state-
elongation and swelling of polydimethylsiloxane networks,” J. Polym. Sci. B-
Polym. Phys, vol. 13, no. 4, pp. 683—702, 1975.

[39] M. C. Boyce and E. M. Arruda, “Swelling and mechanical stretching of elas-
tomeric materials,” Math. Mech. Solids, vol. 6, pp. 641—659, 2001.

184

[40] G. Gee and L. R. G. Treloar, “The interaction between rubber and liquids. 1.
a thermodynimic study of the system rubber-benzene,” Trans. Faraday Soc.,
vol. 38, no. 147—165, 1942.

[41] P. J. Flory, “Thermodynamics of high polymer solutions,” J. Chem. Phys,
vol. 10, pp. 51—61, 1942.

[42] M. L. Huggins, “Some properties of solutions of long-chain compounds,” J. Phys.
Chem, vol. 46, pp. 151-158, 1942.

[43] M. L. Huggins, “Theory of solutions of high polymers,” J. Am. Chem. Soc.,
vol. 64, pp. 1712-1719, 1942.

[44] L. R. G. Treloar, “The equilibrium swelling of cross—linked amorphous polymers,”
Proceedings of the Royal Society of London, vol. A200, pp. 176—183, 1950.

[45] L. R. G. Treloar, “Swelling of a rubber cylinder in torsion: part 1. theory,”
Polymer, vol. 13, pp. 195-202, 1972.

[46] J. J. Shi, K. R. Rajagopal, and A. S. Wineman, “Applications of the theory
of interacting continua to the diffusion of a fluid through a non-linear elastic
media,” Int. J. Eng. Sci, vol. 19, no. 6, pp. 871—889, 1981.

[47] M. V. Gandhi, K. R. Rajagopal, and A. S. Wineman, “Some nonlinear diffusion
problems within the context of the theory of interacting continua,” Int. J. Eng.
Sci, vol. 25, no. 11-12, pp. 1441—1457, 1987.

[48] A. S. Wineman and K. R. Rajagopal, “Shear induced redistribution of fluid
within a uniformly swollen nonlinear elastic cylinder,” Int. J. Eng. Sci, vol. 30,
no. 11, pp. 1583-1595, 1992.

[49] M. V. Gandhi, M. Usman, A. S. Wineman, and K. R. Rajagopal, “Combined
extension and torsion of a swollen cylinder within the context of mixture theory,”
Acta Mech., vol. 79, pp. 81—95, 1989.

[50] S. P. Marra, K. T. Ramesh, and D. A. 8., “Characterization and modeling of
compliant active materials,” J. Mech. Phys. Solids, vol. 51, pp. 1723—1743, 2003.

[51] S. P. Marra, K. T. Ramesh, and D. A. S., “The actuation of a biometic poly(vinyl
alcohol)-poly(acrylic acid) gel,” Phil. Trans. R. Soc. A, vol. 360, pp. 175—198,
2002.

[52] R. W. Ogden, “Large deformation isotropic elasticity: On the correlation of
theory and experiment for incompressible rubberlike solids,” Phil. Trans. R. Soc.
A, vol. 326, pp. 565—584, 1972.

[53] J. Dolbow, E. Fried, and H. Ji, “A numerical strategy for inverstigating the ki-
netic response of stimulus-responsive hydrogels,” Comput. Methods Appl. Mech.
Eng, vol. 194, no. 4447-4480, 2005.

185

[54] J. Dolbow, E. Fried, and H. Ji, “Chemical induced swelling of hydrogels,” J.
Mech. Phys. Solids, vol. 52, no. 51-84, 2004.

[55] J. Dolbow, E. Fried, and H. Ji, “Kinetics of thermally induced swelling of hy-
drogels,” Int. J. Solids Struct, vol. 43, no. 1878-1907, 2006.

[56] S. Ramtani, “Steady diffusion of an ideal fluid through a pre-stressed and rein-
forced hollow conduit subjected to combined finite deformations,” Int. J. Solids
.Struct., vol. 44, no. 14-15, pp. 4819—4829, 2007.

[57] K. Garikipati, E. M. Arruda, and K. Grosh, “A continuum treatment of growth
in biological tissue: the coupling of mass transport and mechanics,” J. Mech.
Phys. Solids, vol. 52, pp. 1595—1625, JULY 2004.

[58] R. M. Bowen, “Theory of mixtures,” Cont. Phys, vol. III, pp. 1—127, 1976.

[59] B. Markert, “A biphasic continuum approach for viscoelastic high-porosity
foams: Comprehensive theory, numerics, and application,” Arch. Comput.
Method Eng, vol. 15, pp. 371—446, DEC 2008..

[60] M. Landervik and R. Larsson, “Modeling of large inelastic deformations of foams
with respect to the point of compaction,” Eur. J. Mech. A-Solids, vol. 27, no. 2,
pp. 234-246, 2008.

[61] P. H. af Segerstada, R. Larssona, and S. Toll, “A constitutive equation for open-
cell cellular solids, including viscoplasticity, damage and deformation induced
anisotropy,” Int. J. Plast, vol. 24, no. 5, pp. 896-914, 2008.

[62] J. Korsawe and G. Starke, “A least-squares mixed finite element method for biots
consolidation problem in porous media,” SIAM J. Numer. Anal, vol. 43, no. 1,
pp. 318—339, 2005.

[63] J. Korsawe, G. Starke, and W. Wang, “Finite element analysis of poro—elastic
consolidation in porous media: Stande and mixed approaches,” Comput. M eth-
ods Appl. Mech. Eng, vol. 195, no. 9—12, pp. 1096—1115, 2006.

[64] B. Markert, “A constitutive approach to 3—d nonlinear fluid flow through finite
deformable porous continua - with application to a high-porosity polyurethane
foam,” Transport Porous Med, vol. 70, no. 3, pp. 427—450, 2007.

[65] W. Ehlers, “Constitutive equations for granular materials in geomechanical con-
text,” Cont. Mech. Environ. Sci. Geophys, vol. 337, pp. 313-402, 1993.

[66] W. Ehlers and G. Eippers, “Finite elastic deformations in liquid-saturated and
empty porous solids,” Transport Porous Med, vol. 34, pp. 179—191, 1999.

[67] A. P. del Palomar and M. Doblare, “On the numerical simulation of the mechan-
ical behaviour of articular cartilage,” Int. J. Numer. Meth. Eng, vol. 67, no. 9,
pp. 1244—1271, 2006.

186

[68] S. Diebels, “A macroscopic description of the quasi-static behavior of granular
materials based on the theory of porous media,” Granul. Matter, vol. 2, pp. 143—
152, JUN 2000.

[69] R. Larsson, M. Wysocki, and S. Toll, “Process-modeling of composites using
two—phase porous media theory,” Euro. J. Mech. A-Solids, vol. 23, pp. 15-36,
JAN-FEB 2004.

[70] B. Markert, B. Monastyrskyy, and W. Ehlers, “Fluid penetration effects in porous
media contact,” Contin. Mech. Therm, vol. 20, pp. 303—315, OCT 2008.

[71] S. Baek and A. R. Srinivasa, “Diffusion of a fluid through an elastic solid un-
dergoing large deformation,” Int. J. Non-Linear Mech., vol. 39, pp. 201—218,
2004.

[72] W. Hong, X. Zhao, J. Zhou, and Z. Suo, “A theory of coupled diffusion and large
deformation in polymeric gels,” J. Mech. Phys. Solids, vol. 56, pp. 1779—1793,
2008.

[73] W. Hong, Z. Liu, and Z. Suo, “Inhomogeneous swelling of a gel in equilibrium
with a solvent and mechanical load,” Int. J. Solids Struct, vol. 46, pp. 3282-
3289, 2009.

[74] R. J. Monroe, “Measuring the mechanical response of swollen hydrogels,” Mas-
ter’s thesis, Michigan State University, 2008.

[75] R. S. Rivlin, “Large elastic deformation of isotropic materials. vi. further results
in the theory of torsion, shear and flexure,” Phil. Trans. R. Soc. Lond. A, vol. 242,
no. 845, pp. 173—195, 1949.

[76] O. H. Varga, “Stress-strain behavior of elastic materials,” Interscience, 1966.

[77] D. M. Haughton and A. Orr, “On the eversion of incompressible elastic cylin-
ders,” Int. J. Non-Linear Mech., vol. 30, pp. 81—95, 1995.

[78] D. M. Haughton and A. Orr, “On the eversion of compressible elastic cylinders,”
Int. J. Solids Struct., vol. 34, pp. 1893—1914, 1997.

[79] D. M. Haughton and A. Orr, “Further results for the eversion of highly com-
pressible elastic cylinders,” Math. Mech. Solids, vol. 1, pp. 355-367, 1996.

[80] Y.-C. Chen and D. M. Haughton, “Existence of exact solutions for the eversion
of elastic cylinders,” J. Elast., vol. 49, pp. 79—88, 1997.

[81] J. G. Simmonds and P. Warne, “Azimuthal shear of compressible or incom-
pressible, non-linearly elastic polar orthotropic tubes of infinite extent,” Int. J.
Non-Linear Mech., vol. 27, pp. 447-464, 1992.

187

 

[82] A. S. Wineman and W. K. Waldron Jr., “Normal stress effects induced during
circular shear of a compressible non-linear elastic cylinder,” Int. J. Non-Linear
Mech., vol. 30, pp. 323—339, 1995.

[83] N. Kikuchi and N. Triantafyllidis, “On a certain class of elastic materials with
non-elliptic energy densities,” Q. Appl. Math, vol. 40, pp. 241—248, 1982.

[84] E. S. Tarleton, J. P. Robinson, and M. Salman, “Solvent-induced swelling of
membranes - measurements and influence in nanofiltration,” J. Membr. Sci,
vol. 280, pp. 442-450, Sept. 2006.

[85] E. S. Tarleton, J. P. Robinson, S. J. Smith, and J. J. W. Na, “New experimen-
tal measurements of solvent induced swelling in nanofiltration membranes,” J.
Membr. Sci, vol. 261, pp. 129—135, SEP 2005.

[86] K. R. Rajagopal, A. S. Wineman, and M. Gandhi, “On boundary conditions for
a certain class of problems in mixture theory,” Int. J. Eng. Sci, vol. 24, no. 8,
pp. 1453—1463, 1986.

[87] T. J. Pence, “Comments on the continuum theory of finitely deformed solid-fluid
mixtures,” unpublished notes.

[88] S. C. Parasad and K. R. Rajagopal, “On the diffusion of fluids through solids
undergoing large deformations,” Math. Mech. Solids, vol. 11, pp. 291—305, 2006.

[89] K. R. Rajagopal and L. Tao, Mechanics of Mixtures. World Scientific, 1996.

[90] M. Buonsanti, R. Fosdick, and G. Royer-Carfagni, “Chemomechanical equilib-
rium of bars,” J. Elast., vol. 84, pp. 167-188, 2006.

[91] H. Demirkoparan and T. J. Pence, “Swelling of an internally pressurized nonlin-
early elastic tube with fiber reinforcing,” Int. J. Solids .Struct., vol. 44, no. 11-12,
pp. 4009—4029, 2007.

[92] G. Y. Qiu and T. J. Pence, “Remarks on the behavior of simple directionally
reinforced incompressible nonlinearly elastic solids,” J. Elast., vol. 49, no. 1,
pp. 1—30, 1997.

188