Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow
(USC Thesis Other)
Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
CHEMICAL AND MECHANICAL DEFORMATION OF POROUS MEDIA
AND MATERIALS DURING ADSORPTION AND FLUID FLOW
by
Sahar Bakhshian
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CHEMICAL ENGINEERING)
May 2018
Copyright 2018 Sahar Bakhshian
Dedication
To Hassan
and my parents.
ii
Acknowledgements
First and foremost, I would like to express my sincere appreciation to my advisor, Professor
Muhammad Sahimi, for his help, his role in advising me, his insights into the research matters,
as well as issues pertaining to and beyond the academic life and, more generally, his support
throughout my Ph.D. work and beyond that, a culmination of which is this dissertation. He is
responsible for not just my scientic achievements, but also my mentorship and education on
various aspects of transport phenomena in chemical engineering, as well as life. His steadfast
quest to understanding and contributing to new avenues in the eld of porous materials have
served as a constant source of inspiration for me.
I would like to thank my Ph.D. committee members, Professors Theodore T. Tsotsis, Katherine
S.-F. Shing, Rajiv Kalia, Aiichiro Nakano, and Felipe de Barros. Their feedback, guidance,
motivating talks and support are truly appreciated.
I would also like to thank Professor Ahmed Elbanna of the University of Illinois at Urbana
Champaign, as well as Professors Kristian Jessen and Theodore T. Tsotsis of USC. It was a
pleasure to work with them. In particular, I thank Professor Tsotsis for providing experimental
data that is used in Chapter 3.
I have a special debt to Dr. Mahnoush Babaei, Professor Sahimi's wife, who helped me during
my serious illness and treatment. I am very grateful for the time she took to talk with me and
my family to help us better understand my treatment options.
iii
Finally, I want to thank my family and friends who helped me and supported me all the way.
In particular, I would like to thank my husband Hassan and my sister Sonia for their unbounded
and sel
ess support.
This work was supported as part of the Center for Geologic Storage of CO
2
, an Energy
Frontier Research Center, funded by the U.S. Department of Energy, Oce of Science, Basic
Energy Sciences, under Award DE-SC0C12504.
iv
Table of Contents
Dedication ii
Acknowledgements iii
List of Tables vii
List of Figures viii
Abstract xiii
Chapter 1: Introduction 1
1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Organization of This Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Chapter 2: Adsorption-Induced Swelling of Clay Minerals 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Numerical Simulation and Model Implementation . . . . . . . . . . . . . . . . . . . 14
2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 Estimating the Model's Parameters . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.2 Sorption Isotherms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.3 Sorption and Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.4 Eect of the Microstructure . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.5 Eect of the Matrix Rigidity . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.6 Eect of Fluid-Solid Interaction . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.5 Possible Implications for Microseismic Events . . . . . . . . . . . . . . . . . . . . . 34
2.6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
Chapter 3: Image-based Modeling of Gas Adsorption and Deformation in Porous
Media with Application to Sandstones 36
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2 Measurement of the Sorption Isotherm . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 Numerical Simulation of Sorption and Deformation . . . . . . . . . . . . . . . . . . 39
3.4 Numerical Calculation of the Eective Permeability . . . . . . . . . . . . . . . . . 42
3.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.5.1 Estimating the Parameters of the Model . . . . . . . . . . . . . . . . . . . . 43
3.5.2 The Permeability of the Core Sample . . . . . . . . . . . . . . . . . . . . . 44
3.5.3 Adsorption and Volumetric Strain Isotherms . . . . . . . . . . . . . . . . . 44
3.5.4 Eect of the Hydrostatic Stress . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.5.5 Adsorption under Geological Conditions . . . . . . . . . . . . . . . . . . . . 51
3.6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
v
Chapter 4: Adsorption and Deformation in Flexible Metal-Organic Frameworks 56
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2 Computer Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3.1 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3.2 Adsorption Isotherms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.3 Phase Fractions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.3.4 Strain and Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Chapter 5: Computer Simulation of Eect of Deformation on the Morphology
and Flow Properties of Granular Porous Media 70
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2 Generation of Random Packings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.3 Equations of Motion and Simulation of Deformation . . . . . . . . . . . . . . . . . 74
5.4 Calculation of the Eective Permeability . . . . . . . . . . . . . . . . . . . . . . . . 78
5.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.5.1 Distributions of the Pores' Size and Length . . . . . . . . . . . . . . . . . . 83
5.5.2 Porosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.5.3 Number of Contacts between Particles . . . . . . . . . . . . . . . . . . . . . 91
5.5.4 Fluid Stress Distribution and Deformation . . . . . . . . . . . . . . . . . . . 92
5.5.5 Eective Permeability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.5.6 Eect of the Friction Coecient . . . . . . . . . . . . . . . . . . . . . . . . 99
5.6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Chapter 6: Eect of Flash Heating on the Frictional Strength Evolution of Gran-
ular Media During Slip Weakening 102
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.2 Theoretical Model and Computer Simulation . . . . . . . . . . . . . . . . . . . . . 104
6.2.1 Discrete-Element Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
6.2.2 Friction Weakening by Flash Heating . . . . . . . . . . . . . . . . . . . . . 107
6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Bibliography 117
vi
List of Tables
2.1 Dependence of the parameters
and
on the microstructure. . . . . . . . . . . . 27
3.1 Mineralogy of the core samples according to their x-ray diraction . . . . . . . . . 38
3.2 Properties of the core samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.1 Dependence of the parameters
and
on temperature for CO
2
. . . . . . . . . . 64
4.2 Dependence of the parameters
and
on temperature for CH
4
. . . . . . . . . . 66
5.1 Numerical values of the exponent m and the prefactor n [see Eq. 5.21] for various
types of packings, particle-size distributions, and type of particles (soft and hard). 88
6.1 Micromechanical and geometrical properties of the packing. . . . . . . . . . . . . . 109
vii
List of Figures
2.1 The microstructure of the system. The white areas represent where adsorption
occurs, while the gray represents the solid skeleton. The triangles are due to nite-
element discretization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 (a) Experimental data
1
(squares) and calculated (dashed line) excess adsorption
isotherm at T
= 0:73. (b) Calculated absolute adsorption isotherm at T
= 0:73
for
=20 and
= 1:38. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3 CO
2
adsorption/desorption isotherms for
= 1:38 at T
= 0:73. Darker regions
represent higher densities of CO
2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Snapshots of CO
2
adsorption along the red isotherm in Figure 3 for P=P
0
= (a)
0.12; (b) 0.14; (c) 0.26; (d) 0.35; (e) =0.56, and (f) 1.0, for
= 0,
= 1:38, and
T
= 0:73. Darker regions represent higher densities of CO
2
. . . . . . . . . . . . . 21
2.5 Snapshots of CO
2
adsorption along the black isotherm in Figure 3 for P=P
0
= (a)
0.12; (b) 0.14; (c) 0.26; (d) 0.35; (e) 0.56, and (f) 1.0, for
=20,
= 1:38, and
T
= 0:73. Darker regions indicate higher densities of CO
2
. . . . . . . . . . . . . . 23
2.6 Dependence of (a) swelling strain, and (b) porosity on the pressure for
= 1:38
and T
= 0:73. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.7 Porosity-dependence of strain for
= 1:38 and T
= 0:73. . . . . . . . . . . . . . . 24
2.8 Variation of strain with CO
2
density for
=20,
= 1:38, and T
= 0:73. . . . 25
2.9 Directional strains for
= 1:38 and T
= 0:73. . . . . . . . . . . . . . . . . . . . . 26
2.10 The microstructures used in the simulations. (a) Triangular; (b) square; (c) Voronoi;
(d) random pore space and (e) parallel platelets. . . . . . . . . . . . . . . . . . . . 28
2.11 Adsorption/desorption isotherms computed with various microstructres and the
parameters shown in Table 1 at temperature T
= 0:73. . . . . . . . . . . . . . . . 29
2.12 Dependence of the strain on the pressure computed with various microstructure
and the parameters shown in Table 1 at temperature T
= 0:73. . . . . . . . . . . 30
viii
2.13 Dependence of the strain on the pressure and the Young's modulusY of the matrix
for
=20,
= 1:38, and T
= 0:73. . . . . . . . . . . . . . . . . . . . . . . . . 31
2.14 Dependence of the adsorption/desorption isotherms on the Young's modulus Y of
the solid, with
=20,
= 1:38, and T
= 0:73. . . . . . . . . . . . . . . . . . . 31
2.15 Adsorption/desorption isotherms for
=20 and T
= 0:73, and their depen-
dence on the
uid-solid energy interaction parameter
. . . . . . . . . . . . . . . . 33
2.16 Dependence of the strain on the pressureP and the
uid-solid interaction parameter
for
=20 and T
= 0:73. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.1 Image of the sandstone sample (left) and its discretization with adaptive tetrahedral
mesh. Green represents that pore space, while blue indicates the solid skeleton. The
image was segmented and discretized using Seg3D 2.4.0 (CIBC, 2016a) and Cleaver
2.2 (CIBC, 2016b), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2 Comparison of the experimental data with the computed adsorption isotherms.
Also shown are the computed desorption isotherms. . . . . . . . . . . . . . . . . . 47
3.3 Computed dependence of the volumetric strains corresponding to the isotherms of
Figure 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4 CO
2
(a) absolute and (b) excess adsorption in the porous sample at 313 K and 350
K. No hydrostatic pressure was applied to the structure. . . . . . . . . . . . . . . . 48
3.5 (a) Volumetric strain and (b) the porosity, corresponding to the isotherms of Figure
3.4, at 313 K and 350 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6 Absolute and excess adsorption, as well as the volumetric strain in the porous
medium under a hydrostatic pressure of 50 MPa and two temperatures. . . . . . . 50
3.7 Same as in Figure 3.6, except that the hydrostatic pressure is 100 MPa. . . . . . . 52
3.8 Strain-dependence of CO
2
excess and absolute adsorption in the porous medium
under three hydrostatic pressures and temperatures of (a,c) 313 K and (b,d) 350 K. 53
3.9 Dependence on the depth of the (maximum) absolute and excess adsorption of CO
2
with a hydrostatic pressure gradient of 15 MPa/km and temperature gradient of
27.3
C/km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1 Schematic of the MIL-53 framework with the cell parameters a,b andc. The view
is along the c axis. (a) Large pore, and (b) narrow pore structures. . . . . . . . . . 59
4.2 Discretized structures of the (a) lp, and (b) the np frameworks, using tetrahedral
cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3 Comparison of computed CO
2
adsorption isotherms with and experimental data. . 61
4.4 Comparison of computed CH
4
adsorption isotherms with the experimental data. . 62
ix
4.5 Pressure-dependence of the fractions of the np and lp phases in the transition zone
during CO
2
adsorption. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.6 Pressure-dependence of the fraction of the np and lp phases in the transition zone
during adsorption of CH
4
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.7 The computed volumetric strain in the MIL-53(Al) during CO
2
adsorption in the
lp and np frameworks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.8 The computed volumetric strain in the MIL-53(Al) during CH
4
adsorption in the
lp and np frameworks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1 The error,jK
N
K
250
j=K
250
in calculating the eective permeabilityK
N
of a grid
of size NNN relative to K for N = 250. . . . . . . . . . . . . . . . . . . . . . 78
5.2 The three particle-size distribution (PSD) used in the for random packings, along
with an example of a random packing before and after compression. . . . . . . . . 79
5.3 The computed pore-size distributions of the dense packings of hard particles with
the three particle-size distributions shown in Fig. 5.2. . . . . . . . . . . . . . . . . 81
5.4 The computed pore-length distributions of the dense packings of hard particles
with the three particle-size distributions shown in Fig. 5.2. . . . . . . . . . . . . . 82
5.5 The computed pore-size distributions for the dense packings of soft particles and
the three particle-size distributions shown in Fig. 5.2. . . . . . . . . . . . . . . . . 84
5.6 The computed pore-length distributions for the dense packings of soft particles and
the three particle-size distributions shown in Fig. 5.2. . . . . . . . . . . . . . . . . 85
5.7 The computed pore-size distributions of random packings of hard particles with an
initial porosity 0.54 and the three particle-size distributions shown in Fig. 5.2. . . 86
5.8 The computed pore-length distributions of the random packings of hard particles
with an initial porosity 0.54 and the three particle-size distributions shown in Fig.
5.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.9 Same as in Fig. 5.7, but for random packings of soft particles. . . . . . . . . . . . . 87
5.10 Same as in Fig. 5.8, but for random packings of soft particles. . . . . . . . . . . . . 87
5.11 Fluctuations in the porosity in the strain (z) and
ow (x) directions. . . . . . . 89
5.12 Dependence of porosity of the dense packings on the external strain. The particle-
size distributions are shown in Fig. 5.2, and the porosities have been normalized
by their initial values before deformation. . . . . . . . . . . . . . . . . . . . . . . . 89
5.13 Same as in Fig. 5.12, but for loose packing with an initial porosity of 0.54. The
insets show the logarithmic plots of the same. . . . . . . . . . . . . . . . . . . . . . 90
x
5.14 Same as in Fig. 5.12, but for the regular packings. The insets show the logarithmic
plots of the same. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.15 Dependence of the average number of contacts C on the external strain in random
packings of soft particles and the particle-size distributions shown in Fig. 5.2. . . . 91
5.16 The distributions of the computed normalized stresses, exerted by the
uid in the
pore space of the dense random packings of soft particles. The curves represent t
of the numerical data to the log-normal distribution, Eq. 5.23. . . . . . . . . . . . 92
5.17 The distribution of the computed normalized stresses, exerted by the
uid in the
pore space of a loose packing of soft particles with an initial porosity of 0.54. The
curve show the t of the numerical data to Eq. 5.23. . . . . . . . . . . . . . . . . . 93
5.18 The distribution of the computed normalized stresses, exerted by the
uid in the
pore space of regular packings of soft equal-size particles. (a) Simple-cubic (SC)
packing at an external strain of 0.2, and (b) face-centered cubic (FCC) packing at
an external strain of 0.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.19 Dependence of the permeability of dense random packings on the external strain.
The particle-size distributions are those shown in Fig. 5.2. The permeabilities are
normalized by their values before deformation. . . . . . . . . . . . . . . . . . . . . 95
5.20 Dependence of the permeability of dense random packings on the porosity. The
particle-size distributions are those shown in Fig. 5.2. The permeabilities are
normalized by their values before deformation. . . . . . . . . . . . . . . . . . . . . 96
5.21 Collapse of the permeability-porosity data of Fig. 5.20 onto a universal curve. . . . 97
5.22 Same as in Fig. 5.19, but for random loose packings with an initial porosity of 0.54. 97
5.23 Dependence of the permeability of the random loose packings on the porosity. The
initial porosity was 0.54. The particle-size distributions are those shown in Fig.
5.2. The permeabilities are normalized by their values before deformation. . . . . . 98
5.24 Collapse of the permeability-porosity data of Fig. 5.23 onto a universal curve. . . . 99
5.25 Dependence on the external strain of the permeability of the simple-cubic (SC) and
face-centered cubic (FCC) packings of particles. The permeabilities are normalized
by their values before deformation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.26 Dependence on the friction coecient of the permeability and porosity of a ran-
dom packing of soft particles with an initial porosity of 0.54. The particle-size
distribution is PSD3, shown in Fig. 5.2. . . . . . . . . . . . . . . . . . . . . . . . . 100
6.1 Random packing of spheres used in the shear simulation. The upper-layer particles
move with velocityV
s
in thex direction, while the lower-level particles are held xed.107
6.2 Shear stress-strain for constant and variable friction factors with shear velocity =
1 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
xi
6.3 Shear stress-strain for constant and variable friction factors with shear velocity =
0.5 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.4 Shear stress-strain for constant and variable friction factors with shear velocity =
0.2 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.5 Shear stress-strain for constant and variable friction factors with shear velocity =
0.1 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.6 Shear stress-strain for constant and variable friction factors with shear velocity =
0.08 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.7 Shear stress-strain for constant and variable friction factors with shear velocity =
0.05 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.8 Dependence of steady shear stress on shear velocity for weakening and no-weakening.114
6.9 Distribution of weakened friction factor for various shear velocities and at times
(all in s)(a) 2.5; (b) 500; (c) 1000, and (d) 2150. . . . . . . . . . . . . . . . . . . . 115
xii
Abstract
Theoretical analysis and numerical modeling of morphology of porous materials and their
ow,
transport, and deformation have been widely studied for decades, due to the wide number of
applications of such materials in petroleum, geotechnical, and chemical engineering, materials
science, and biomechanics. An important class of porous materials is those that are deformed
by, for example, a large external pressure or shear, or by
uid adsorption that causes swelling
and/or shrinkage of the porous medium. The research presented in this thesis deals with distinct
modes of deformation in porous media and its eect on the
uid
ow and transport properties
of the pore space. One type is concerned with deformation that arises as a result of injecting a
gas or liquid into the pore space, and the resulting induced deformation that is the result of the
interaction between the solid matrix of the medium and its
uid content.
In the second type we study deformation of a porous medium as a result of applying an external
strain to the medium, and investigate the resulting changes in the morphology of the pore space,
namely, its pore-size distribution, pore-length distribution, pore connectivity and porosity, as well
as its eective permeability.
We rst study deformation of porous media induced by sorption of
uids, and in particular by
carbon dioxide. The phenomena are signicant, not only to the capacity of geological formations
for storing CO
2
, but also to the eect that they have on the mechanical properties of the formations
and microseismic events that they might possibly trigger. To study the phenomena, we formulate
the problem by energy consideration in which the Hamiltonian (total energy) of a porous medium
and its
uid content is represented as the sum of the elastic energy of the system, as well as
xiii
the energies associated with the interactions between the adsorbates and between them and the
medium's solid matrix. The gas phase is described by the density-functional theory, while the
solid matrix can be either linearly or nonlinearly elastic. The model provides predictions for
the sorption isotherms, the dependence of the strain on the bulk pressure, and the change in
the porosity. When the model was applied to sorption of CO
2
in clay particles, all the reported
experimental features of the phenomena were reproduced.
We also use the same model for adsorption of CO
2
and its induced deformation in a sandstone
formation over wide temperature and pressure ranges. Finite-element computations are used to
compute CO
2
adsorption isotherms and the induced strain in the formation. All the computations
are carried out with three-dimensional image of a core sample from Mt. Simon sandstone, the
target formation for a pilot CO
2
sequestration project being conducted by Illinois State Geological
Survey. The computed CO
2
sorption isotherm at 195 K is in excellent agreement with the available
experimental data. The model is capable of predicting gas adsorption in porous formations at
high pressures and temperatures. Thus, it is used to study the eect of hydrostatic pressure on
adsorption and deformation of the porous formation under various conditions. We nd that the
eect of the conning pressure is more prominent at higher temperatures. Also computed is the
depth-dependence of the capacity of the formation for CO
2
adsorption, along with the induced
volumetric strain.
Next, we utilize the same statistical mechanical model of adsorption-induced deformation of
porous media in order to study the same phenomena in
exible metal-organic frameworks (MOFs).
MOFs have recently attracted considerable attention as new nanoporous materials with applica-
tions to separation, catalysis, gas capture and storage, and drug delivery. A subclass of MOFs,
the MIL-53(Al or Cr) family, exhibits a complex structural phase transition, often referred to as
the breathing transition, which occurs when gas molecules are adsorbed in its pore space. In this
phenomenon, the materials' morphology oscillates between two distinct phases, usually referred
to as the narrow-pore (np) and large-pore (lp) phases. We use the model to study adsorption of
xiv
CO
2
and CH
4
in MIL-53(Al) samples. The results demonstrate that all the reported experimental
features of the phenomenon are produced by the model. In particular, consistent with experimen-
tal data, the model predicts that in the np phase, the material experiences contraction, whereas
in lp phase it undergoes swelling. In the pressure range in the transition zone, the coexistence of
the two phases is demonstrated.
We also report on the results of extensive computer simulation of the eect of mechanical
deformation on the morphology of a granular porous medium and its
uid
ow properties. The
porous medium is represented by packings of spherical particles. Both random and regular, as
well as dense and loose packings are used. A quasistatic model based on Hertz's contact theory is
used to model the mechanical deformation of the packings, while the evolution of the permeability
with the deformation is computed by the lattice-Boltzmann approach. The evolution of the pore-
size and pore-length distributions, the porosity, the particles contacts, the permeability, and the
distribution of the stresses that the
uid exerts in the pore space are all studied in detail. The
distribution of the pores lengths, the porosity, and the particles connectivity change strongly with
the application of an external strain to the porous media, whereas the pore-size distribution is not
aected as strongly. The permeability of the porous media strongly reduces even when the applied
strain is small. When the permeabilities and porosities of the random packings are normalized
with respect to their predeformation values, they all collapse onto a single curve, independent of
the particle-size distribution. The porosity reduces as a power law with the external strain. The
uid stresses in the pore space follow roughly a log-normal distribution, both before and after
deformation.
In the last chapter of this thesis we model shear strength loss at frictional asperity contacts,
induced by
ash heating or
ash weakening in a granular fault gouge, using the discrete-element
method. The shear test simulations were carried out at slip rates ranging from 0.05 to 1 m/s.
This study is an eort to demonstrate the velocity-dependent frictional weakening in granular
media. The response of local friction coecients to variation of sliding rate has been studied. We
xv
studied the magnitude of shear-stress reduction and local friction factor variation at various shear
velocities. The results indicate that the local frictional coecients decay during the weakening.
For higher shear velocities, the local friction coecients between the particles take on values lower
than the static friction coecient. It is also found that the frictional stress decreases by applying
the
ash heating model in the shear test simulations. At velocities higher than 0.1 m/s, the eect
of frictional weakening is discernible in the frictional strength (shear-stress strain behavior) of the
medium. The results demonstrate that the eect of frictional weakening is more prominent at
higher sliding velocities.
xvi
Chapter 1
Introduction
1.1 Background and Motivation
Theoretical analysis and numerical modeling of the morphology of porous materials and the phe-
nomena that occur there have been widely studied for decades, due to their wide applications
in petroleum engineering, geotechnical engineering, chemical engineering, materials science, and
biomechanics. An important class of porous materials includes those that are deformed by, for
example, a large external pressure or shear, or by
uid adsorption or desorption that cause de-
formation - swelling or contraction. Important examples of deformable porous materials include
consolidating clay and shale formations (Sahimi, 2011; Bakhshian and Sahimi, 2017), biopolymers
such as cellulose, hydrogels, and cell membranes (Iritani et al., 2006; Karada, 2010; Salimi et al.,
2010; Savoji and Pourjavadi, 2006), paper (Ghassemzadeh and Sahimi, 2004; Masoudi and Pillai,
2010), snow (Wu et al., 2005; Carbone et al., 2010), fabric ltration (Sahimi and Imdakm, 1991),
aggregates formed by
occulation (Meakin, 1998), foams (Koehler et al., 2000; Hutzler et al., 2004;
Pitois et al., 2009; Lorenceau et al., 2009), soft biomaterial such as articular cartilage (Lai and
Mow, 1980), drug delivery systems (Leonowicz and Oleszak, 2009), metal-organic frameworks for
gas storage that will be studied later in this thesis (Mason et al., 2014; He et al., 2014; Li et al.,
2012), and soil and rock at high depths subjected to large external pressure or shear deformation
1
(Sahimi, 2011; Biot, 1941, 1956; Iliev et al., 2008; Bakhshian and Sahimi, 2016; Ferdowsi et al.,
2013; Scuderi et al., 2015).
Deformation of porous materials due to their liquid content, or by storage of gases in their pore
space, is a widespread phenomenon. Print papers swell when they come into contact with water,
or even by the coating
uid that makes high quality surface (Ghassemzadeh et al., 2001). In
biomedical science, the interactions between ion sand charged proteoglycans govern the deforma-
tion of cartilaginous soft-hydrated tissues. Swelling of polymers and biopolymers have numerous
applications, ranging from size-exclusion chromatography, to ecient materials for drug delivery,
gel electrophoresis during ltration processes, contact lenses, and food stu.
An important class of porous media to which deformation is important is geomaterials. For
example, porous formations that contain charged molecules and brine may swell if penetrated by
fresh water. Moreover, swelling of smectite clays by water adsorption has been studied for decades.
Clay swelling is of critical importance to geotechnical and geoenvironmental elds, as it may cause
damage to the foundations of buildings, giving rise to cracks and displacements of footings. These
are irreversible phenomena with considerable economic implications. Bentonitic-based compacted
clays have low hydraulic conductivity, plasticity (Benson et al., 1994), swelling and adsorptive
capacity for contaminants, which is why they are used as sealing materials to prevent leaking of
contaminants to the environment. Storage of radioactive nuclear wastes in geological formations
is perhaps the only viable solution to problem of safeguarding such materials, and compacted
clays may be used for sealing the storage area. Yet another example is provided by gas-bearing
shales that represent unconventional sources of energy that have attracted wide attention. Their
swelling is a major negative factor during drilling, as it can cause signicant problems for wellbore
stability.
Fluid adsorption causes deformation of porous media, which occurs at molecular scale through
the interactions between solid and
uid molecules (?). The deformation of porous structure
induced by adsorption of molecules or ion has been considerably studied in recent years. It has
2
been shown that the structural change of solid matrix is an inseparable and signicant part of the
adsorption process. Since adsorption-induced deformation in most of nanoporous materials has
been observed to reach the order of 1%, it is neglected in most of the cases (Dang et al., 2017).
One example for this phenomenon is deformation of coal seams and shale formations upon gas
adsorption during storage of CO
2
and enhanced recovery of coalbed methane (Zang and Wang,
2017; Mazumder and Wolf, 2008; Bakhshian and Sahimi, 2017). Although the deformation is
very small, the associated stresses may mount up to high values and induce large external forces
on the formations. Moreover, adsorption-induced swelling of underground formations has been
identied as one of the main factors that aects porosity, permeability, diusivity and surface area
of the pore space (Siriwardane et al., 2009; Pan et al., 2010). Another example is metal-organic
frameworks (MOFs), a relatively new class of microporous materials used mainly for gas storage.
These materials can be subjected to structural transformation and deformation of large amplitude
in response to
uid adsorption (Suh et al., 2012; Jeroy et al., 2008).
There have been many experimental studies aimed at understanding the response of the porous
materials during adsorption (Cui et al., 2007; Chen et al., 2012; Serre et al., 2007; Lyu et al.,
2015; van Bergen et al., 2009). In addition, various theoretical approaches have been developed
to predict the stress induced by adsorption and its eect on solid matrix deformation in porous
materials (Bakhshian and Sahimi, 2017; Vandamme et al., 2010; Brochard et al., 2012). Despite
such eorts, a clear picture of how the phenomenon occurs at the pore scale has not emerged.
A part of this thesis deals with simulation of CO
2
adsorption and its induced deformation
in geological formations, such as shale and sandstone reservoirs. We describe a comprehensive
theoretical framework to describe adsorption-induced deformation in a porous medium with a
complex pore space geometry. Since adsorption in a porous medium with disordered geometry
poses dicult challenges to its modeling, most of the previous studies considered a regular pore
geometry, such as slit or cylindrical pores, or a collection of independent simple pores. Such
models do not take into account the complex interconnectivity of the pores, a characteristic of all
3
porous media. Thus, the proposed model is an attempt to study
uid adsorption and the resulting
induced deformation in disordered porous materials.
On the other hand, deformation of porous media due to the pore pressure and external me-
chanical forces is also an important problem in geology. Reservoir rocks and their pore structures
are constantly under stresses arising from various triggers, including the pore
uid or the conn-
ing pressure from the surrounding rock formations (Neuzil, 2003; Mair et al., 2002; Saksala, 2016;
Ougier-Simonin and Zhu, 2013). The eld stresses may promote rock deformation, compaction
or failure in the subsurface medium. Depletion of hydrocarbon reservoirs during production and
injection of high volume of CO
2
during its sequestration in geological formations are some of the
problems that induce stress on the formations and cause mechanical deformation (Zoback and
Zinke, 2002; Vilarrasa et al., 2016). Such deformations may induce compaction of the formations,
decrease in their permeability and porosity, give rise to wellbore stability issues, sand production,
subsidence and fault reactivation (Fossen et al., 2007; Alam et al., 2014; Dautriat et al., 2009).
Since unconsolidated rocks have lower shear and tensile strength and higher compressibility than
the consolidated ones, they are more susceptible to the eects of mechanical deformation.
Due to the importance of this phenomenon in oil and gas industry, it is imperative to un-
derstand how microstructural,
ow, and transport properties of a formation vary in response to
mechanical deformation. Numerous experimental and numerical studies have been conducted to
understand the mechanical response of rock formations while they are under compression. Numer-
ical modeling has advantages over experimental studies due to their capabilities for investigating
the behavior of the system under various conditions that are dicult to experiment in laboratory
setting. Thus, a part of this thesis focuses on the computer simulation of the eect of mechanical
deformation on the morphology of a porous medium and its
uid
ow properties.
4
1.2 Organization of This Thesis
Emission of greenhouse gases, such as CO
2
produced essentially by burning fossil fuels, to the
atmosphere is the main contributor to the global warming. In recent years, sequestration of CO
2
in
geological formations is considered as one option to mitigate the undesirable environmental impact
of CO
2
emission to the atmosphere. A signicant fraction of CO
2
can be stored in geological
formations by adsorption, which is accompanied by medium deformation - swelling and shrinkage.
Thus, knowledge of adsorption behavior of CO
2
and its induced swelling is important in conducting
feasibility studies of CO
2
storage in subsurface reservoirs. Despite tremendous experimental
studies by previous researchers over the past decades, a modeling approach that can simulate the
coupling of adsorption and its induced deformation in porous structures is still scarce. This subject
has motivated us to propose a comprehensive model by coupling the lattice gas density functional
theory (DFT) with the equation for elastic deformation and nite-element (FE) computations in
order to study adsorption-induced swelling in a disordered porous medium. We have used the
model to study sorption of CO
2
in clay minerals and its eect on their swelling (Chapter 2).
In the current theoretical and computational methods that are utilized to study CO
2
sorption
in rocks, a simple pore space model is utilized that has practically no relation with the actual
shapes and sizes of the pores. For example, molecular simulation studies of CO
2
sorption in
rock have utilized a slit pore, the space between two
at and parallel surfaces. To address this
shortcoming, we have used our theoretical framework to simulate sorption and its associated
deformation in three-dimensional image of a core sample from Mt. Simon sandstone formation
(Chapter 3).
As pointed out earlier, carbon capture is a way for tackling the problem of global warming
and emission of greenhouse gases, such as CO
2
. Among the carbon capture technologies, solid
sorption oers a promising method for ecient and selective sorption of CO
2
with low energy cost,
compared with common and dominant amine-based absorption systems. While sorbents such as
5
zeolite and activated carbon have been widely used for CO
2
capture applications, signicant energy
required for their regeneration and low CO
2
loading are of disadvantages of these materials.
An intriguing class of porous materials that are known as metal-organic frameworks (MOFs)
are synthesized and studied as alternative sorbents for the purpose of gas storage and separation.
Some classes of such materials represent noticeable stimuli-responsive behavior, such as abrupt
structural deformation upon gas adsorption. Up to now, there have been very limited number
of computational studies to investigate the behavior of such materials. To better understand
gas adsorption and its induced structural transformation in MOFs, we study these phenomena
through our adsorption-deformation model (Chapter 4).
Pore pressure drop during oil production and injection of high volume of CO
2
during its
sequestration in geological formations induce large stresses exerted on the formations, which have
a signicant eect on the microstructural and petrophysical properties of reservoir rocks. Since
such eects are very important for oil and gas industry, this subject has motivated us to study the
eect of mechanical deformation, induced by
uid
ow, on the morphology of a granular porous
medium and its
uid
ow properties through extensive computer simulations (Chapter 5).
Another phenomenon that causes deformation of rock formations is attributed to the fault
zone during sliding. Fault zones contain granular cores, known as fault gouge. Earthquakes
occur when a fault gouge weakens with increasing slip rates. Thus, fault gouge plays an essential
role in its shear strength and the slip dynamics of earthquake. Understanding the mechanism
of shear weakening in a granular medium is an underlying challenge in the earthquake physics.
Hence, inspired by
ash weakening process during earthquakes, we have used the discrete-element
method to model shear strength loss at frictional asperity contacts, induced by
ash heating or
ash weakening in a granular gouge (Chapter 6).
6
Chapter 2
Adsorption-Induced Swelling of Clay Minerals
2.1 Introduction
As described in Chapter 1, one problem in which swelling of porous media is important is seques-
tration of carbon dioxide in geological formations, such as depleted oil and gas reservoirs, and
unminable coal seams (Michels et al., 2015; de Jong et al., 2014; Hemmen et al., 2012; Loring
et al., 2012; Tsotsis et al., 2004; Yang et al., 2011; Billemont et al., 2013; McDonald et al., 2015;
Cho et al., 2015; Dadwhal et al., 2008), as well as in metal-organic frameworks and other types
of adsorbents that have been fabricated for capturing CO
2
(Sumida et al., 2012). The main char-
acteristics of the geological formations for long-term CO
2
geo-sequestration include high storage
capacity and eective sealing of the caprocks that are composed frequently of shale and, thus, clay
minerals. The low permeability of the clay formations reduces the possibility of upward migration
of CO
2
in the formations and its penetration into the neighboring natural resources (Shukla et al.,
2010)
Clay minerals are typically composed of smectite, kaolinite, chlorite and illite (Loring et al.,
2014; Gasparik et al., 2014), among which smectites have the highest sorption capacity due to
their high internal surface area (Ji et al., 2012). Busch et al. (2008) studied CO
2
sorption in a
shale sample and various clay minerals in order to investigate the storage capacity and sealing
7
integrity of such formations. They reported that clay minerals, especially montmorillonite, have
signicant sorption and storage capacity for CO
2
.
More generally, numerous experimental studies have been carried out to understand adsorption
of CO
2
in porous media. For example, a study by Ravikovitch et al. (2005) on adsorption of CO
2
in soil samples containing various organic matters at 273 K and low pressures indicated that CO
2
tends to adsorb onto the domain of soil that is rich in organic matter. In another experimental
study of CO
2
storage in organic-rich shales, it was found that even though such formations are
nearly impermeable and have low porosity, they are capable of storing a remarkable amount of
CO
2
, which was attributed to the nely-distributed organic nanopores (Kang et al., 2011). Many
careful and precise experiments, as well as theoretical studies have also been reported on the eect
of CO
2
sorption in coal (Vandamme et al., 2010; Hol et al., 2011, 2012; Hol and Spiers, 2012).
Theoretical and computational studies of
uid sorption in porous media have also been carried
out, including lattice density-functional theory (DFT) (Pan and Connell, 2007; Ottiger et al.,
2008; Woo et al., 2001), quenched solid DFT (Yang et al., 2011; Gor and A. V. Neimark, 2011),
nonlocal DFT (Ravikovitch and Neimark, 2006), thermodynamic models (Grosman and Ortega,
2008; Kulasinski et al., 2015b; Liu et al., 2015), molecular dynamics (MD) (Ghassemzadeh et al.,
2000; Xu et al., 2001; Ghassemzadeh and Sahimi, 2004; Kim et al., 2007, 2005; Liu and Bhatia,
2013; Rao and Leng, 2016) and Monte Carlo simulations (Valiullin et al., 2006).
Adsorption of gases in a porous matrix leads to its swelling, which aects the porosity, perme-
ability, diusivity, and surface area of the pore space (Perera et al., 2011; Ozdemir and Schroeder,
2009), and has been studied both experimentally and theoretically (Ravikovitch and Neimark,
2006; Grosman and Ortega, 2008; Balzer et al., 2011). For example, a combined experimental and
modeling study of adsorption of CO
2
and N
2
in coal samples has been reported, which demon-
strated that volumetric strain caused by gas adsorption leads to reduction of the permeability
(Pini et al., 2009). Another experimental study of coal swelling in response to CO
2
adsorption
by Anggara et al. (2014) demonstrated that coal tends to swell anisotropically. The maximum
8
volumetric swelling ranged from 1.60 to 2.35 percent. Recent MD simulation of diusion and
adsorption of CO
2
in Montmorillonite structures as the model of clay minerals indicate that the
basal spacing of the material changes. Pan and Connell (2007) developed a theoretical framework
based on energy consideration to describe volumetric changes in coal due to gas adsorption. Their
model was capable of predicting changes in the porosity and permeability of coalbeds as a result
of swelling or shrinkage by gas adsorption. Vandamme et al. (2010) developed a model capable
of calculating the macroscopic strains in response to adsorption in a porous medium. Using a
thermodynamic framework, they developed the poromechanics of the process by considering the
surface energy. The surface stress was modied by adsorption, a phenomenon that has also been
studied by molecular dynamics simulation. A thermodynamic framework was also developed by
Hol et al. (2012) that modeled adsorption of CO
2
by unconned and conned coal. They com-
bined their model with poroelasticity and studied the competition between adsorption-induced
strain and poroelastic strain in terms of the applied external stress,
uid pressure and concentra-
tion. In a study by Espinoza et al. (2013) a poromechanical model for adsorption and strain was
developed. Their model predicted that strain anisotropy is essentially due to anisotropy in the
mechanical properties of porous solid.
Several studies of swelling induced by CO
2
sorption have been reported in coal samples (Yang
et al., 2011; Kim et al., 2011; Kowalczyk et al., 2010), as well as in clay minerals and shales
(Giesting et al., 2012; Rother et al., 2013; Kadoura et al., 2016; Rao and Leng, 2016). One such
study was reported by Heller and Zoback (2014) who carried out experiments on the adsorption
capacity of CO
2
and CH
4
of shale samples from four reservoirs, as well as pure illite and kaolinite.
Their study indicated that gas adsorption causes a volumetric swelling strain on the order of 10
5
to 10
3
. In a more recent study swelling of shales as a result of CO
2
adsorption was investigated.
All the samples studied exhibited anisotropic strain caused by CO
2
adsorption (Lu et al., 2016).
In terms of deforming rock - but not the mechanism by which the deformation is achieved -
sorption-induced deformation of porous materials is somewhat similar to subsidence of oil and gas
9
reservoirs (Doornhof et al., 2006). The mechanism of the two are, of course, dierent. Sorption-
induced deformation is due to the chemical potential of the sorbing species, whereas subsidence is
due gas pressure reduction. Flow of oil or gas in some reservoirs causes deformation, compaction
and subsidence in the reservoir. Thus, to model such phenomena one must couple the elastic or
elasto-plastic equations governing the deformation to the equations of hydrodynamics governing
uid
ow. In a similar spirit, any reasonable model of sorption-induced swelling of a porous
medium must couple the elastic deformation of the solid to the thermodynamic state of the
uid
in the pore space, the chemical potential of the sorbing species, and the interactions between the
solid and the gas (Chin and Nagel, 2004; Guyer and Kim, 2015; Vandamme et al., 2010; Hol et al.,
2012; Espinoza et al., 2013).
The goal of this chapter is to present and use one such model by coupling the lattice gas DFT
with the equation for elastic deformation and nite-element (FE) computations in order to study
adsorption-induced swelling in a disordered porous medium, and in particular clay particles. We
use the model to study sorption of CO
2
in clay minerals and its eect on their swelling.
The rest of this chapter is organized as follows. In the next section we describe the details of
the model for calculating the adsorption isotherms and the deformation strains. Next, the details
of the numerical simulations are presented. The model is then applied to CO
2
sorption in clay and
the resulting swelling and change in its porosity. The results are then presented and discussed.
In particular, we study the eect of the model's parameters on sorption, strain and porosity. The
extension of the model to changing the mechanical properties of porous media, as well as possible
triggering microseismic events that have recently been reported, are discussed. We summarize the
key conclusions in the last section.
10
2.2 The Model
Let us rst state clearly our assumptions. We assume, (i) a quasi-static system, and (ii) a linearly-
elastic solid matrix for the porous media (see also below). The coupling between the solid matrix
and its deformation and the
uid, and the interactions between the two are formulated through
the total energy E of the system (the porous medium plus its
uid content). Thus,
E =E
f
+E
s
+E
fs
; (2.1)
where E
f
, E
s
, and E
fs
denote, respectively, the energy due to the
uid, the solid, and the
interactions between the two. Various expressions can be used for each of the three terms on the
right side of Eq. (2.1). One of the simplest of such expressions for E
f
is given by a mean-eld
DFT (Pitois et al., 2009; Kierlik et al., 1998, 2001),
E
f
=
1
X
i
[
i
ln
i
+ (1
i
) ln(1
i
)]
1
2
ff
X
i
X
Z
i
i+Z
X
i
i
; (2.2)
where = (k
B
T )
1
, withk
B
being the Boltzmann's constant, and T the temperature. The basis
of Eq. (2.2) is a lattice-gas model in which each site i is occupied by only one species, and only
nearest-neighbor interactions are allowed. Thus,Z is the coordination number of the structure or
the lattice used. In Eq. (2.2) is the external chemical potential,
i
is the average
uid density at
sitei, and
ff
is the
uid-
uid interaction parameter. Each sitei is characterized by two indicator
functionsI
i
andI
i
, such thatI
i
= 0 and 1 and (1I
i
) = 0 and 1 represent respectively, the
uid
and the matrix indicator functions, so that
i
=hI
i
I
i
i.
Since we assume that the solid matrix is linearly elastic, the expression for E
s
is straightfor-
ward:
E
s
=
1
2
X
i
(K
i
"
2
i
+ 2
i
"
i
); (2.3)
11
where K
i
is the elastic stiness, and "
i
and
i
are the strain and stress at i, respectively. For
linearly-elastic solids the relation between the strain and displacement is known and used here.
But, the assumption of elastic linearity can also be relaxed, in which case Eq. (2.3) must be
replaced with a constitutive equation for an elastically-nonlinear solid (Sahimi, 2003; Guyer and
Johnson, 2009). It remains to specify the energy E
fs
associated with the
uid-solid interactions.
Clearly, E
fs
consists of two terms:
E
fs
=E
w
+E
ce
: (2.4)
Here,E
w
is associated with the \wetting energy" that is present at any
uid-solid interface, while
E
ce
is the change in the elastic energy as a result of the presence of the
uid. Thus, we write
(Guyer and Kim, 2015),
E
w
=
1
2
fs
X
i
i
; (2.5)
E
ce
=
1
2
X
i
[
i
i
(r u)
i
]; (2.6)
where u is the displacement vector,
fs
is a
uid-solid interaction parameter related to the wetting
properties of the surface, while
i
is a parameter associated with the change in the strain energy
and deformation of the solid as a result of the
uid's presence. The sum is over all pairs of
the lattice sites. The porous medium and its content are in
uid exchange with a reservoir at
temperature T and chemical potential .
The physical meanings of the interaction parameters
fs
and
ff
are clear. They can be
estimated by a mean-eld theory (see below), or by direct MD computations. But, at rst sight,
it might not be clear how to estimate the parameter
i
by MD simulations. Although in this
study we treat
i
as an adjustable parameter, it can in fact be computed by MD calculations
(eorts in this direction are underway), the procedure for which is as follows. Let L be the length
of the pore or void in which the
uid resides. To estimate
i
one needs some sort of representation
12
of the pore, either as a collection of the atoms properly bound together to represent the solid, as
in standard MD simulations, or a continuum representation of the pore. The mechanical state
of the structure must be precisely and properly specied. Then, the
uid is introduced into the
pore or void that is in equilibrium with the chemical potential reservoir (bulk state). The pore
walls have an attractive potential that draws the
uid particles toward them. Thus, in principle,
a well-dened thermodynamic state with energy E(T;V;) exists, whereV is the volume. As the
chemical potential is increased, the
uid's congurations change and try to contract/expand the
pore length L. One monitors and computes the variations in L as a function of . Then, as Eq.
(2.6) indicates,
i
= dL=d(), or equivalently,
i
= dL=dN, where N is the number of particles
on pore wall.
It would perhaps be useful to compare our approach with the previous works that also at-
tempted to couple sorption with deformation. Prominent among such models is the work of Hol
and co-workers. Hol et al. (2012) consider the Gibbs energy of the solid and the
uid phase and
add them together to obtain an expression for the total energy of the system. By dierentiating
the total free energy of the system and considering the condition of equality of chemical potential
of adsorbed CO
2
and the free CO
2
in the equilibrium state, they derived an equation that repre-
sents the equilibrium concentration of the adsorbed CO
2
. Then, by changing the
uid pressure at
a constant temperature, the adsorbed concentration of CO
2
was computed at various pressures.
Finally, the volumetric strain induced by the adsorption was calculated by using the adsorbed
concentration of CO
2
and the molar volume of the adsorbed CO
2
.
To obtain the equilibrium state of the system, we minimize the total energy E with respect to
the density and strain (displacement) to derive the governing equations for the local density and
the displacement eld. This is described next.
13
2.3 Numerical Simulation and Model Implementation
For a given realization of a porous medium and a specic bulk chemical potential that is in
equilibrium with the pore space and its
uid content, the set of equations resulting from minimizing
the total energy with respect to u and
i
must be solved numerically. The minimization yields
the following equation for the
uid density
i
,
i
=
(
1 + exp
[ +
fs
+
ff
X
Z
i+Z
i
(r u)
i
]
!)
1
: (2.7)
For a given chemical potential we rst provide guesses for
i
and u and substitute them in
the right side of Eq. (2.7), which yields a new approximate solution for
i
, which enables us to
discretize the second governing equation by the FE method using triangular elements and linear
basis functions. The second equation (after minimizing the total energy) for the displacement u
l
of vertex l in triangle i is given by,
k
e
u
l
=
1
2
e
ni
i
i
; (2.8)
where the e
ni
are the three unit vectors from the center of triangular element i to its vertices n.
Note that the right side of Eq. (2.8) is simply the force acting on the vertices of the triangular
elements. The domain of the computations is a clay particle, including its pore and solid phases;
see Fig. 2.1.
To proceed with computing the displacements, we rst calculate the stiness matrix of each
FE triangular element. If a triangle's vertices are at (x
1
,y
1
), (x
2
,y
2
), and (x
3
,y
3
), then its stiness
matrix k
e
is given by
k
e
=
Z
S
B
T
DBdS =SB
T
DB; (2.9)
14
Figure 2.1: The microstructure of the system. The white areas represent where adsorption occurs,
while the gray represents the solid skeleton. The triangles are due to nite-element discretization.
withS =
1
2
J =
1
2
(x
13
y
23
x
23
y
13
) being the triangular element's area, D the elasticity matrix,
B the strain-displacement matrix, and T denoting the transpose operation. Here, x
ij
=x
i
x
j
and similarly for y
ij
. Then,
B =
1
J
0
B
B
B
B
B
B
@
y
23
0 y
31
0 y
12
0
0 x
32
0 x
13
0 x
21
x
32
y
23
x
13
y
31
x
21
y
12
1
C
C
C
C
C
C
A
; (2.10)
and
D =
Y
(1 +)(1 2)
0
B
B
B
B
B
B
@
1 0
1 0
0 0
1
2
[1 2]
1
C
C
C
C
C
C
A
; (2.11)
whereY is the Young's modulus (bulk modulus for the
uid phase), and the Poisson's ratio. We
then set up a global stiness matrix K by compiling in it the stiness matrices of the triangular
elements, which gives rise to the following system of linear equations describing the displacements
of elemental vertices:
X
l
K
nl
u
l
=
1
2
X
i
e
ni
i
i
: (2.12)
15
The set 2.12 is solved by the conjugate-gradient method to obtain a new approximate solution for
the displacements. We then substitute the new approximations in the right side of Eq. (2.7) to ob-
tain the next approximate solution for
i
, return to Eq. (2.12) to compute the new displacements,
and continue the iteration until
i
and the displacements no longer change between the iterations.
Periodic boundary conditions were used for computing the density eld. As for the displacement
eld, we used a free-displacement boundary condition. This was meant to mimic the conditions
under which the experimental data with which we compare the model's predictions were measured.
In the experiments (Michels et al., 2015) the sample was in contact with a bath of CO
2
, unlike
a reservoir condition under compaction with an external stress applied to the external surface of
the particle. The corresponding porosity of the pore space and the pressure-dependence of the
strain are then computed. Note that each FE triangular element is either lled by the
uid or
represents the solid matrix. To begin the iteration process for solving the set of equations (2.7)
and (2.12), a uniform density prole was used as the initial approximation. In order to compute
the adsorption isotherms, the chemical potential was increased incrementally from a very large
negative value to the bulk chemical potential. In the present work we consider 2D disordered
porous media, as the experimental data that we used to determine some of the parameters of the
model were for a thin, essentially 2D medium. But, it would pose no particular diculty to use
a 3D model. The number of the triangles representing the solid matrix was selected so that the
porosity of the model pore space matches that of the porous medium with which we compare the
results.
2.4 Results
We have used the model to study a variety of cases. In what follows we present and discuss the
results.
16
2.4.1 Estimating the Model's Parameters
The model is general in that, given expressions for the various types of energies described in
the last section, it can be applied to adsorption-induced swelling in any porous material. Our
main interest in this chapter is, however, the adsorption of CO
2
in clay particles. As usual,
the sorption isotherms are obtained by computing the overall density of CO
2
in the pore space,
CO2
= N
1
P
i
i
as a function of the pressure P at a temperature below the bulk critical
temperature, with N being the total number of sites in the pore space. In addition, we compute
the variations of the strain and the porosity during sorption. The excess adsorption is calculated
by,
ex
=
1
N
X
i
(
i
b
); (2.13)
where
b
is the density of the bulk
uid. The properties of the bulk
uid are given by (Monson,
2012),
p
N =Nk
B
T ln(1
b
)
1
2
ZN
ff
2
b
; (2.14)
b
=k
B
T ln
b
1
b
Z
ff
b
; (2.15)
where
p
= Pv
s
is the lattice equivalent of the bulk pressure, and v
s
is the physical volume
equivalent to a lattice site. Equations (2.14) and (2.15) are for the lattice-gas model.
The interaction parameters
ff
and
fs
, the strain energy parameter
i
, and the maximum
pore-lling capacity that is required for converting the physical units of the experimental data for
sorption to lattice gas units must be determined. According to the mean-eld theory, the bulk
critical temperature is given by
T
c
=
Z
ff
4k
B
; (2.16)
17
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
P/P
0
ρ
ex
Model
Experiment
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P/P
0
ρ
CO
2
(b)
(a)
Figure 2.2: (a) Experimental data
1
(squares) and calculated (dashed line) excess adsorption
isotherm at T
= 0:73. (b) Calculated absolute adsorption isotherm at T
= 0:73 for
=20
and
= 1:38.
and, thus, the
uid-
uid interaction parameter is given by,
ff
=k
B
= 4T
c
=Z. With CO
2
being
the adsorbate we obtain,
ff
=k
B
405:67 K, if we use T
c
304:25 and Z = 3 for CO
2
. The
maximum pore-lling capacity is n
max
v
p
, where n
max
is the maximum
uid density, and v
p
is
the pore volume. We rst estimate the range of variations of n
max
. The minimum value of
n
max
is estimated by the inverse of the van der Waals co-volume, which for CO
2
is 1.03 gr/cm
3
.
Its maximum is estimated as
p
2=
3
, the density of a close packing of spheres, where is the
molecular diameter. Therefore, the upper limit ofn
max
for CO
2
is (Woo et al., 2001) 1.74 gr/cm
3
.
Thus,n
max
varies by a factor less than 2. In fact, the theoretical sorption isotherm is not sensitive
to n
max
. Our preliminary simulations indicated that the lower limit of n
max
provides accurate
results and, therefore, we xed its value at the lower theoretical limit. The pore volumev
p
is xed
by the experimental values of the density and porosity of the experimental samples, which for the
experimental data (Michels et al., 2015) that we used were 2.8 gr/cm
3
and 0.4, respectively. The
remaining parameters,
fs
and
i
, were determined by tting the calculated adsorption isotherm
to the data.
18
The experimental data that we used were for excess CO
2
adsorption by a smectite clay (syn-
thetic Li-
uorohectorite) sample at room temperature and pressures varying from 0 to 45 bars.
The sample had a volume of 20 mm
3
, with a thickness of 1 mm. The elastic modulus and the
Poisson's ratio of the material were taken to be, respectively, 2 GPa and 0.2 (although higher and
lower values of the modulus have also been reported (Bathija, 2009)). Since the porous samples
used in the experiments were thin, we used a 2D model whose microstructure, a random distri-
bution of polygons representing the solid matrix in a uniform background, is shown in Figure 2.1.
The eect of the microstructure is studied later. The model's dimensions were selected such that
they match those of the experimental sample and its porosity of 0.4 and, as pointed out earlier,
periodic boundary conditions were used in the computations for the density eld. The calcula-
tions with the FE method involved 2476 triangular element. We checked that there is no nite-size
eect by carrying out simulations with a larger system. The initial computations were carried
out in order to determine the optimal values of the remaining two parameters. We obtained,
=
fs
=
ff
1:38 and
=
i
=
ff
20 for the dimensionless experimental temperature,
T
= k
B
T=
ff
0:73. Unless specied otherwise, these are the values that we utilized in every
case.
2.4.2 Sorption Isotherms
Figure 2.2(a) presents a comparison between the computed isotherm and the experimental data,
and the agreement is excellent. Both the experiments and the computed isotherm indicate that
with increasing pressure, the excess adsorbed CO
2
increases smoothly to a maximum value, and
then decreases at higher pressures. Figure 2.2(b) represents the corresponding absolute sorption
isotherm, which reaches a plateau and remains unchanged at higher pressures.
Before studying various aspects of the model and its predictions, it is worth looking at the
eect of the parameter
i
on the isotherms. Recall that Eq. (2.6) represents the change in the
elastic energy as a result of the presence of the adsorbed
uid with
i
0. Thus,
i
provides a
19
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P/P
0
ρ
CO
2
−20
0
λ
*
Figure 2.3: CO
2
adsorption/desorption isotherms for
= 1:38 at T
= 0:73. Darker regions
represent higher densities of CO
2
.
link between the
uid adsorption and the solid deformation. When
i
= 0, deformation of the
solid is not aected by adsorption. Since a t of the theoretical predictions to the data yielded,
20, we also studied the case,
= 0, holding all other parameters xed.
In all the cases that we presented hereafter and the corresponding gures, the solid curves
present the state of the dependent variable during adsorption, while the dashed curve indicate the
same during desorption. Figure 2.3 presents the adsorption/desorption isotherms for two values,
= 0 and20. The existence of a plateau in the sorption isotherm suggests, of course, that
there is an optimal pressure for maximum storage of CO
2
. If
= 0, i.e. when sorption does
not aect solid swelling, sorption hysteresis is very weak, almost indiscernible. As the absolute
value of
increases, coupling deformation of the solid with adsorption, the
uid forces the solid
to work harder in order to accommodate more
uid. Thus, hysteresis becomes stronger and more
visible, because deformation changes the pore structure that is responsible for the hysteresis. The
shape of the isotherms for the two cases are, however, identical, implying that the type of sorption
20
Figure 2.4: Snapshots of CO
2
adsorption along the red isotherm in Figure 3 for P=P
0
= (a) 0.12;
(b) 0.14; (c) 0.26; (d) 0.35; (e) =0.56, and (f) 1.0, for
= 0,
= 1:38, and T
= 0:73. Darker
regions represent higher densities of CO
2
.
21
isotherms and the pore-lling mechanism are not aected by the swelling. The deformation only
assists more adsorption. The maximum amount of adsorbed
uid is higher when swelling and
sorption are coupled. Figure 2.4 presents the snapshots of the evolution of CO
2
density prole
in the porous medium with increasing pressure for
= 0, i.e. when sorption and deformation
are decoupled. Darker areas correspond to higher amounts of adsorbed CO
2
. This should be
compared with Figure 2.5 that presents the same evolution, but for the case,
=20. It is clear
that the FE triangular elements that are neighbors to the solid skeleton of the porous medium
are occupied by CO
2
before the other sites, which once again demonstrates the mutual eects of
sorption and deformation.
2.4.3 Sorption and Deformation
Figure 2.6 presents the dependence of the total volumetric strain and the porosity of the system
on the pressure. Due to the rigidity of the matrix the porosity increases by only about 3.4 percent
over its original value, and even less when sorption is decoupled from deformation (
= 0). This
can also be understood by appealing to the variations of the total strain on the pressure, also
shown in Figure 2.6. If we plot the strain versus the porosity, we obtain the results shown in
Figure 2.7, indicating a linear relation between the two, which for small strains is expected and
indicates that swelling is caused by gas adsorption (Heller and Zoback, 2014; Sahimi, 2011). The
fact that the slope of the line is greater than one indicates that the matrix itself has swelled,
contributing to the total strain (Kulasinski et al., 2015a).
Figure 2.8 presents the swelling strain as a function of adsorbed amount of CO
2
, indicating a
linear relationship between the two. Although we do not have such data for the clay samples that
we used for estimating the parameters of the model, we note that a linear relation between the
strain and the adsorbed amount is in agreement with the experimental data for CO
2
adsorption in
shales at low pressures (Lu et al., 2016). Day et al. (2008) also reported that coal swelling and the
amount of CO
2
adsorbed into it are linearly correlated for pressures less than 8 MPa. Liu et al.
22
Figure 2.5: Snapshots of CO
2
adsorption along the black isotherm in Figure 3 for P=P
0
= (a)
0.12; (b) 0.14; (c) 0.26; (d) 0.35; (e) 0.56, and (f) 1.0, for
=20,
= 1:38, and T
= 0:73.
Darker regions indicate higher densities of CO
2
.
23
0 0.2 0.4 0.6 0.8 1
0.402
0.404
0.406
0.408
0.41
0.412
0.414
0.416
P/P
0
Porosity
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
P/P
0
Strain (%)
−20
0
−20
0
(a) (b)
λ
*
λ
*
Figure 2.6: Dependence of (a) swelling strain, and (b) porosity on the pressure for
= 1:38 and
T
= 0:73.
0.4 0.405 0.41 0.415 0.42
0
0.005
0.01
0.015
0.02
0.025
Porosity
Strain
0.4 0.405 0.41 0.415 0.42
0
0.005
0.01
0.015
0.02
0.025
Porosity
Strain
(a)
λ
*
=0
(b)
λ
*
=−20
Figure 2.7: Porosity-dependence of strain for
= 1:38 and T
= 0:73.
24
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
ρ
CO
2
Strain (%)
Figure 2.8: Variation of strain with CO
2
density for
=20,
= 1:38, and T
= 0:73.
(2015) studied the volumetric strain in unconned bituminous coal samples, exposed to adsorbing
CO
2
at pressures up to 100 MPa and a temperature of 40
C. Their study also conrmed a linear
relation between adsorption-induced swelling strain and the amount of adsorbed CO
2
.
To better illuminate the eect of sorption on the strain and deformation of porous media,
we computed the longitudinal and transverse strains,
x
and
y
, respectively. For example, the
longitudinal strain
x
, the change in the length of the matrix relative to its original length L
0
in
the x direction (horizontal direction in Figure 2.1), is dened by
x
=
x
L
0
: (2.17)
A similar equation is used for
y
. Figure 2.9 presents the two strains and their dependence on
pressure and the parameter
i
. We rst note that strain changes anisotropically, the reason for
which is twofold. One is the deviation of the microstructure of the system from elastically-isotropic
structures, such as hexagonal and triangular lattices. The second reason is due to the dierences
25
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
P/P
0
Directional Strain (%)
ε
x
ε
y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
P/P
0
Directional Strain (%)
ε
x
ε
y
(b) λ
*
=−20 (a) λ
*
=0
Figure 2.9: Directional strains for
= 1:38 and T
= 0:73.
between the dierent elastic properties of the matrix and the
uid content of the pores. If we
dene a strain anisotropy by,A =
y
=
x
, then its value at the maximum pressure P =P
0
is 1.4
for
= 0, and 1.36 if
< 0. Although we do not have any data forA for the experimental clay
system with which we compare our results, our computedA is similar to the experimental data
for CO
2
adsorption in shale samples for which the average anisotropy was reported to be 1.41
(Lu et al., 2016). In addition, Figure 2.9 illustrates hysteresis in the strain as the porous medium
rst adsorbs CO
2
and then desorbs it. As with the sorption isotherms, the hysteresis is stronger
when sorption is fully coupled with deformation of the solid matrix, since it is the structure of
the material that gives rise to hysteresis, and swelling changes the microstructure.
2.4.4 Eect of the Microstructure
It is of interest to study the eect of the microstructure of a porous medium on the sorption
isotherms of CO
2
and the resulting swelling of the solid matrix. In particular, it is well-known
that adsorption is not sensitive to pore interconnectivity, whereas desorption is, which is also
why in many cases there is hysteresis between the two (Sahimi, 2011). Thus, we carried out a
series of simulations with various microstructures. They included the cases in which the solid
26
matrix was represented by square and triangular tiles, as well as the Voronoi polygons. We also
considered the case in which the particle is represented by parallel pores sandwiched between
platelets representing the solid matrix. All the microstructures are shown in Figure 2.10.
We rst used the experimental data (Michels et al., 2015) to t the two parameters
and
for the Voronoi, triangular, square and parallel platelet microstructures. The results are
listed in Table 1, and have a few notable features. (i) The parameters for the Voronoi and
triangular structures are practically identical. The pore connectivity of the triangular network
is six, as is the average connectivity of the pores in the Voronoi microstructure. Moreover, the
values of the parameters for the random microstructure, which has a mean connectivity very close
to 4, and those of the square microstructure are also identical. Thus, as long as the average
connectivity of a disordered microstructure is the same as, or very close to, that of a regular
one, the parameters of the model are essentially the same, an important result that was pointed
out long ago (Jerauld et al., 1984a,b) in the context of conduction in disordered microstructures.
(ii) The energy interaction parameter
varies by greatly as the microstructure changes. Table
1 indicates that
for the square microstructure is about 30 percent smaller than that of the
parallel platelet microstructure, which has the highest value. (iii) The parameter
exhibits
stronger dependence on the microstructure than does
, but it still varies by a factor less than
2. Such observations indicate the possibility of developing mean-eld theories for estimating
and
, or computing them by MD simulations.
Table 2.1: Dependence of the parameters
and
on the microstructure.
Random Square Triangular Voronoi Parallel Platelets
1.38 1.38 1.65 1.65 1.98
-20 -20 -12 -12 -20
27
Figure 2.10: The microstructures used in the simulations. (a) Triangular; (b) square; (c) Voronoi;
(d) random pore space and (e) parallel platelets.
Figure 2.11 compares the sorption isotherms computed with the various microstructures. The
dierence between the various isotherms is mainly in the desorption isotherms. In particular, the
hysteresis is most pronounced with the random and square microstructures that also happen to
have the lowest connectivity. This is in line with what was pointed out earlier in that desorption
and, hence, its hysteresis with adsorption are sensitive to the connectivity of a porous material.
Since a pore cannot release the gas during desorption, unless it is connected to other pores that can
do the same at the same pressure, hysteresis should be stronger in less connected microstructures.
Note that the parallel platelets model exhibits no hysteresis, which is expected because the pores
do not communicate with each other.
Figure 2.12 compares the pressure-dependence of the overall strain computed with the various
microstructures. Similar to the sorption isotherms, the hysteresis between the strains build-up
during adsorpion and its release during desorption is most pronounced with the random and
square microstructures. Moreover, Figure 2.12 indicates that the strain accumulation and release
is far more sensitive to the connectivity than are the sorption isotherms. In agreement with Figure
2.11, the strain \isotherms" for the pairs of the microstructures that have similar connectivities
are also close. Once again, the parallel platelets model does not exhibit any hysteresis.
Sorption isotherms have been used in the past in order to estimate the mean connectivity
of porous materials (Seaton, 1991; Liu et al., 1992). Since the pressure-dependence of strain
manifests strong sensitivity to the connectivity, the implication is that one may be able to gain
28
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P/P
0
ρ
CO
2
Random
Voronoi
Parallel Platelets
Triangle
Square
Figure 2.11: Adsorption/desorption isotherms computed with various microstructres and the
parameters shown in Table 1 at temperature T
= 0:73.
more information about the connectivity through mechanical measurements as a porous medium
undergoes deformation during sorption.
2.4.5 Eect of the Matrix Rigidity
We also studied the eect of the rigidity of the solid matrix on deformation and sorption. Clearly,
the more compliant the matrix, the more amenable it is to swelling during sorption. At the same
time, the elastic moduli of clay materials vary over a relatively wide range, depending on their
type and content (impurity). Thus, we carried out a series of simulations in which we assumed
that the Young's modulus of the solid matrix was 0.15 GPa, which was reported for a clay-sand
mixture and is more than ten times less than what we used in most of our calculations (Vallejo
and Lobo-Guerrero, 2005). Also studied was the case in which the modulus was 10 GPa, ve
times larger than the base case.
29
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
P/P
0
Strain(%)
Square
Random
Parallel Platelets
Voronoi
Triangle
Figure 2.12: Dependence of the strain on the pressure computed with various microstructure and
the parameters shown in Table 1 at temperature T
= 0:73.
The resulting pressure-dependence of the strain, its variations with the Young's modulus, and
the comparison with the base case are shown in Figure 2.13. The signicant feature of Figure 2.13
is that the hysteresis between the overall strain during adsorption and desorption decreases as the
matrix becomes more rigid. This is understandable because when the matrix is more compliant,
it is deformed more strongly. The deformation changes the structure of the matrix, which is the
root cause of the hysteresis. Figure 2.14 presents the corresponding sorption isotherms. Similar
to the strain, the hysteresis in the sorption isotherms is much stronger in the more compliant
material. This is of course due to the strong deformation of the matrix that changes the pore
space signicantly, and in particular the pore connectivity, and hysteresis between adsorption and
desorption is caused by the pores' interconnectivity. Moreover, the change in the rigidity also
changes the amount of sorbed gas at low and intermediate pressures. As the matrix becomes
more compliant, it requires less work to accommodate more
uids (see the earlier discussions),
resulting in higher amounts of sorption at low and intermediate pressures.
30
0 0.2 0.4 0.6 0.8 1
0
1
2
3
4
5
6
P/P
0
Strain (%)
0.15
2
10
Y (GPa)
Figure 2.13: Dependence of the strain on the pressure and the Young's modulus Y of the matrix
for
=20,
= 1:38, and T
= 0:73.
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P/P
0
ρ
CO
2
0.15
2
10
Y (GPa)
Figure 2.14: Dependence of the adsorption/desorption isotherms on the Young's modulus Y of
the solid, with
=20,
= 1:38, and T
= 0:73.
31
Interestingly, if we t the two parameters of the model, namely,
and
, to the experimental
data (Michels et al., 2015), assuming the lower or the higher values of the Young's modulus than
the base case, we obtain, respectively,
1:45 and 1.38, which are practically identical with
1:38 that we obtained with the base case. This supports our contention that the interaction
energy parameter between the
uid and the solid is essentially independent of the rigidity of solid
matrix.
2.4.6 Eect of Fluid-Solid Interaction
The interaction between a
uid and the solid is clearly one of the most important factors in
sorption in porous media. The same is true for sorption-induced deformation. Thus, it is of
interest to study the eect of the interaction parameter
fs
on the sorption isotherms and the
change in the strain. The results for three values of
=
fs
=
ff
are presented in Figure 2.15.
The rest of the model's parameters are the same as those of the baseline model, so that a direct
comparison can be made. As
decreases, the hysteresis loop becomes wider, because it is linked
with the thermodynamic metastability that represents a system trapped in a local free-energy
minimum (Gelb et al., 1999). Qualitatively similar behavior and hysteresis loops are seen in the
dependence of the overall strain on the pressure; see Figure 2.16.
The eect of temperature on sorption is well understood, and needs no explanation. As T
increases, the amount of adsorbed gas decreases, implying that we should expect the overall, as
well as the directional strains to also decrease because the swelling is due to sorption. Indeed, our
model does predict this aspect correctly, which is also in agreement with the experimental studies
on CO
2
sorption in shale and coal (Yang et al., 2011; Lu et al., 2016).
32
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P/P
0
ρ
CO
2
0.80
1.38
2.00
ε
*
Figure 2.15: Adsorption/desorption isotherms for
=20 andT
= 0:73, and their dependence
on the
uid-solid energy interaction parameter
.
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
P/P
0
Strain (%)
0.80
1.38
2.00
ε
*
Figure 2.16: Dependence of the strain on the pressureP and the
uid-solid interaction parameter
for
=20 and T
= 0:73.
33
2.5 Possible Implications for Microseismic Events
In addition to their implications for storage of CO
2
, deformation and swelling of porous media,
such as rock, have two important consequences. One is the change that they induce in the
mechanical properties of porous media by swelling. In a sense, the swelling causes the matrix
to soften up, becoming more compliant (although some have disputed this; see, for example,
(Vandamme et al., 2010; Hol et al., 2011, 2012; Espinoza et al., 2013). The second eect is
associated with microseismic events (small earthquakes) triggered by deformation of porous media
and rock that have been reported recently (Rothert and Shapiro, 2003; Maxwell et al., 2008; Tafti
et al., 2013; Lee et al., 2016).
The observed microseismic events may be triggered by one or both of the following mechanisms
(Busch et al., 2008). (i) The deformation and swelling of the matrix activate dormant fractures
and cleavages that already reside in the matrix. Swelling might, of course, also produce an increase
in the normal stress components that can just as well tend to close fractures and stabilize faults
in a formation undergoing CO
2
uptake; see de Jong et al. (2014). The activation of the fractures
and cleavages triggers a microseismic event because deformation changes the stress eld in rock,
and may increase slip along the activated fractures. (ii) On the other hand, swelling softens up
the matrix, making it more amenable to nucleation of new fractures. This is particularly so when
CO
2
is injected into the pore space, which builds up the pore pressure. Once new fractures have
been nucleated, they may propagate and eventually trigger microseismic events.
It has been reported that such events do not usually occur at or very close to the injection
point, but rather at a distance from it. This supports the hypothesis that microseismic events
may be triggered by one or both of the aforementioned mechanisms, because either the dormant
fractures closest to the injection point are rst activated, leading to the activation of other dormant
fractures at a distance from the injection, or new fractures must rst be nucleated and propagated
for some distance, changing the stress eld and leading to a microseismic event. Both indicate
34
that sorption (and
ow that is not considered in this study) of CO
2
may have profound eect
on the structure of the geological formation in which it is to be stored. Thus, one must develop
a comprehensive model that takes into account such eects, in order to better understand what
happens to a geological formation when CO
2
is stored in it. Work in these directions is in progress.
2.6 Summary and Conclusions
This chapter presented a model for deformation and swelling of porous media induced by their
interaction with a
uid, and in particular its sorption. The model is based on the total energy
(Hamiltonian) of the system that takes into account the interaction between the solid matrix of a
porous medium and the adsorbing
uid, as well as the elastic energy stored in the material, and
the interaction energy between the adsorbates themselves. The gas phase is described by a density
functional theory. The model is capable of providing predictions for the sorption isotherms, the
dependence of the strain on the pressure, as well as other properties of interest. The change
in the strain is anisotropic due to the deformation of the solid matrix and the deviation of its
underlying structure from two-dimensional elastically-isotropic materials, as well as the dierence
between the stiness of the matrix and the bulk modulus of the
uid phase. As the adsorbed gas is
desorbed, the strain is also released and, similar to the sorption isotherms, develops hysteresis with
the strain increase during adsorption. This opens up the possibility of mechanical measurements
during sorption experiments for deducing insight into the structure of the porous medium. When
the model was applied to sorption of CO
2
in clay particles, all the reported experimental features
of the phenomena were reproduced by the model. The implications of the results for CO
2
storage
in geological formations were also discussed, and future research directions were suggested.
35
Chapter 3
Image-based Modeling of Gas Adsorption and Deformation
in Porous Media with Application to Sandstones
3.1 Introduction
As described in Chapters 1 and 2, among the many methods that have been proposed for mitigat-
ing the eect of emission of CO
2
, its sequestration in geological formations is considered (Dowell
et al., 2017; Cai et al., 2017; Dashtian et al., 2017; Bakhshian and Sahimi, 2016; Shaer, 2010;
Merey and Sinayuc, 2016; Busch et al., 2008; Liu and Wilcox, 2012) to be a promising way of
addressing the problem of global warming. The geological formations of interest may include
depleted oil and gas reservoirs, deep saline aquifers, and unminable coal seams (Lucier et al.,
2006; Eliebid et al., 2017; Kim et al., 2017). In addition to its sequestration in depleted forma-
tions, the injected CO
2
also provides sustained pressure that enhances the production of CH
4
,
which competes with CO
2
(Heller and Zoback, 2014; Gensterblum et al., 2013) for sorption on
the internal surface of the reservoirs. Better understanding of the phenomenon entails having two
main ingredients. One is the knowledge concerning the adsorption capacity of rocks. The second
ingredient is a model of the formation itself.
Extensive studies of gas adsorption in geological formations have been carried out by both
experiments (Gasparik et al., 2014; Chareonsuppanimit et al., 2012; Loring et al., 2014; Michels
36
et al., 2015; Tsotsis et al., 2004; Day et al., 2008) and various theoretical and computational
methods (Bakhshian and Sahimi, 2017; Zhang et al., 2016; Rao and Leng, 2016; Guo et al., 2016;
Botan et al., 2010; Yang et al., 2011; Psarras et al., 2017; Ghassemzadeh et al., 2000). Grand-
canonical Monte Carlo simulation (Rao and Leng, 2016; Botan et al., 2010) has been the main
computational approach for studying the sorption phenomenon, while such theoretical models as
simplied local-density (Chareonsuppanimit et al., 2012) and the Ono-Kondo models (Guo et al.,
2016; Sudibandriyo et al., 2011; Qajar et al., 2016), as well as the density-functional theory (Yang
et al., 2011; Balzer et al., 2016) constitute the main theoretical methods. In addition, molecular
dynamics simulations (Yang and Zhang, 2005) have provided useful insight into the mechanism
of gas sorption in microporous media. But, despite providing a basis for better understanding
of gas sorption in porous media, almost all the theoretical and computational methods that
have been utilized so far in order to study CO
2
sorption in actual rock suer from two main
shortcomings. One is that models of the pore space that have been utilized have been too simple,
having practically no relation with the actual shapes and sizes of the pores. For example, molecular
simulation studies of CO
2
sorption in rock have utilized a slit pore (Botan et al., 2010; Kadoura
et al., 2016; Yang et al., 2015), the space between two
at and parallel surfaces. The second
shortcoming of the previous studies is the fact that experiments (Heller and Zoback, 2014; de Jong
et al., 2014; Chen et al., 2015; Rother et al., 2013) and limited computational studies (Yang
et al., 2011; Vandamme et al., 2010) have indicated that sorption of CO
2
in rock may swell the
formation, hence changing the structure of the pore space as adsorption proceeds. Adsorption-
induced swelling of porous formations is a main factor that contributes to the change in the
porosity, permeability, diusivity, and surface area of a pore space (Perera et al., 2011).
The main goal of this chapter is to address the aforementioned shortcomings. We use the
statistical mechanical model of gas sorption and the associated deformation in a porous medium
that we developed in Chapter 2 to compute the sorption isotherms and the deformation that it
induces in a core sample from a sandstone reservoir - the Mount (Mt.) Simon formation. The
37
Table 3.1: Mineralogy of the core samples according to their x-ray diraction
Whole-rock mineralogy %
Quartz 73
K-feldspar 12
Plagioclase 7
Ankerite/Fe-dolomite 1
Fluorapatite 1
Pyrite 1
Clay 4
Clay %
Illite-smectite 68
Illite 31
Chlorite 1
Kaolinite 1
formation is generally referred to the basal Cambrian formation in the upper Mississippi Valley and
southern Great Lakes areas. Numerical simulation of the theoretical model is carried out using a
three-dimensional (3D) image of a core sample. We then compare the results with the experimental
sorption isotherms and permeability, measured with the same sample porous medium whose image
we utilize in the simulations, and measured in the research group of Professor Theodore T. Tsotsis
of USC.
3.2 Measurement of the Sorption Isotherm
As mentioned earlier, gas adsorption in two cylindrical core samples from Mt. Simon formation
were measured at USC. The samples are both horizontal plugs from a larger section that had been
extracted along a well in the formation. More specically, sample 1 is from a depth of 6979.57
ft, whereas sample 2 is from a depth of 6985.85 ft. Their dimensions and weights are listed in
Table 3.2. The mineralogy of the parent core sample from which the two cores were extracted is
presented in Table 3.1. Prior to the sorption measurements, the samples' porosities were measured
via the helium expansion method. They are virtually identical, and are listed in Table 3.2.
38
Table 3.2: Properties of the core samples
Core Depth (ft) L (cm) D (cm) W (g) Porosity (%)
1 6979.6 4.88 2.51 46.05 24.8
2 6985.6 5.16 2.47 49.04 24.9
Before sorption experiments began, the samples were dried at 110
C under vacuum (1 mm
Hg) for 24 hours. The N
2
and CO
2
adsorption isotherms were measured using a Micromeritics
ASAP 2010 Instrument. The N
2
sorption test was carried out at the nominal temperature of 77.3
K, which was maintained by immersing the specialty-made sample-holder containing the core
sample in a liquid-nitrogen bath. Note, however, that the true adsorption temperature varies a
bit during the experiment. The instrument, however, measures the true experimental temperature
and calculates the corresponding N
2
saturation pressure. The CO
2
adsorption measurements were
carried out at 195 K, which was maintained by immersing the sample-holder containing the core
in a bath containing a mixture of dry ice and acetone. To generate the adsorption isotherm, the
pressure in the sorption cell was increased in a step-wise manner, and the adsorbed gas volume
was measured and recorded. The data will be compared with the predictions of the model (see
Results and Discussions section).
3.3 Numerical Simulation of Sorption and Deformation
Assuming that the sorption-induced swelling is slow and quasi-static, we use the statistical me-
chanical model, whose details were described in Chapter 2, Section 2.2. We then used the nite-
element method to solve the governing equations. Minimizing the total energy with respect to
i
and u
i
yields two equations, one for the local
uid density
i
, and a second one for the dis-
placement eld of the solid. The two equations are coupled and must be solved numerically. The
governing equations were given in Chapter 2, Section 2.3 (see Eqs. 2.7 and 2.12). We explain
brie
y how to solve the equations using nite-element method for a three-dimensional problem, as
39
this is an important aspect of the problem that was not considered in Chapter 2. The elemental
stiness matrix, k
e
, is given by
k
e
=
Z
V
U
T
CUdV ; (3.1)
where C and U are, respectively, the elastic coecient and the strain-displacement matrices, and
T denotes the transpose operation. The integration is over the volume V of the element in the
computational grid. C is given by,
C =
0
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
@
0
+ 2
0
0
0
0 0 0
0
0
+ 2
0
0
0 0 0
0
0
0
+ 2
0
0 0 0
0 0 0
0
0 0
0 0 0 0
0
0
0 0 0 0 0
0
1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A
; (3.2)
where
0
and
0
are the Lam e constants that are related to the shear modulusG and the Poisson's
ratio ,
0
=
2G
1 2
; (3.3)
0
=G: (3.4)
The 3D computational grid that we utilize consists of tetrahedron grid blocks; see Fig. 3.1.
We rst introduce three normalized coordinates (;;) (f
1
;f
2
;f
3
) that are related to the
Cartesian coordinates (x;y;z) by, x = x
1
+
P
4
i=2
(x
i
x
1
)f
i1
, y = y
1
+
P
4
i=2
(y
i
y
1
)f
i1
,
and z = z
1
+
P
4
i=2
(z
i
z
1
)f
i1
, with (x
i
;y
i
;z
i
) (i = 1 4) being the coordinates of vertices
of the tetrahedra. We then introduce a shape vector, N = (N
1
;N
2
;N
3
;N
4
)
T
whose components
40
N
i
(;;) are, N
1
= 1
P
j = 2
4
f
j1
, and N
j
=f
j1
with j = 2 4. In Eq. 3.1 the matrix U
is a 3 4 matrix given by
U = DN; (3.5)
where D is the dierentiation operator matrix, so that the entries U
ij
of U are given by, U
ij
=
@N
i
=@f
j
. A global stiness matrix K is then constructed by compiling in it the stiness matrices
of the tetrahedral elements. Doing so generates a system of linear equations that describe the
displacements of elemental vertices, which is given by Eq. 2.12.
We compute the sorption isotherms of CO
2
and N
2
and the deformation that they induce in
the 3D core sample from Mt. Simon sandstone, whose actual 3D image of the sample, shown in
Fig. 3.1, we utilize in the simulations. Due to the limitations on the resolution of the image,
the porosity of the image is 10 percent. Thus, based on the correlation function of the pore sizes
(Knackstedt et al., 1998, 2001) we generated pores with sizes below the resolution of the image
and distributed them in the image, so that its total porosity matched that of the core sample.
To carry out the computations, we must rst discretize the 3D image in order to generate the
computational grid. To do so, one must rst determine the size of the representative elementary
volume (REV) for the sample, i.e. the minimum image size such that its properties will not change
if larger images and computational grids are used. In a previous study (Tahmasebi et al., 2016) in
which
ow of CO
2
and water in the same image was simulated, it was determined that an image
with 300 300 300 voxels represents adequately the REV. Thus, we use the same image size in
the present study.
The computational grid must be resolved enough to accurately represent the image's morphol-
ogy, but it must also be such that the required computations do not take prohibitively long times.
We discretized the domain with an adaptive computational grid with tetrahedral elements, shown
in Fig. 3.1. After some preliminary simulations in which various grid resolutions were utilized
in order to identify the most accurate grid resolution with aordable computational time, a grid
41
with 373,607 tetrahedral elements was determined to be accurate enough, as the porosity and
computed adsorption isotherm of CO
2
did not change signicantly if grids with higher resolutions
were used.
3.4 Numerical Calculation of the Eective Permeability
We utilized the lattice-Boltzmann (LB) method with a single-relaxation time to simulate
uid
ow in the pore space of the image. The eective permeability was computed using Darcy's
law. Using the standard bounce-back rule, the no-slip boundary condition was imposed on the
solid boundaries. A constant pressure gradient was applied in the
ow direction. The D3Q19
propagation scheme (Succi, 2001) was used in LB method in order to evolve the
ow system.
The LB model is second-order accurate, with its compressibility error related to the square of
the Mach number, Ma. In order to simulate an incompressible
uid, Ma is kept below 0.1. A unit
relaxation parameter
R
= 1 of the LB method was used, because it leads to accurate simulation
of
uid
ow (Succi, 2001). The simulations were carried out in the creeping
ow regime with
a Reynolds numbers, Re < 1. Some preliminary simulations were carried out to determine an
accurate grid resolution and size. We determined that a grid of size 300300300 yields accurate
estimate of the permeability with reasonable computation time.
3.5 Results and Discussion
In what follows we present and discuss the results, the rst of which is estimation of the parameters
of the model, which we carried out by a method dierent from what was described in Chapter 2.
42
3.5.1 Estimating the Parameters of the Model
The Hamiltonian (total energy) of the system contains three parameters, namely, the
uid-
uid
interaction energy parameter
ff
, the
uid-solid interaction energy parameter
fs
, and the cou-
pling coecient
i
that links the deformation to sorption. One way of estimating
ff
is through
the mean-eld theory that relates the bulk critical temperature T
c
to
ff
by, T
c
= Z
ff
=(4k
B
),
where k
B
is the Boltzmann's constant. As is well-known, however, mean-eld theories usually
over- or underestimate various properties of
uid-solid systems. A more accurate way of esti-
mating
ff
and
fs
is based on the energy interaction parameters that are used in molecular
dynamics (MD) simulations. As Table 3.1 indicates, x-ray diraction of the core samples of Mt.
Simon shows that quartz - silica - constitutes about three-fourth of the sandstone's composition.
In fact, the quartz content of several core samples from the same location in Mt. Simon ranges
from 73 percent to over 90 percent. Therefore, we use the Lennard-Jones (LJ) energy parameters
for CO
2
and CO
2
-silica (Ravikovitch et al., 2001; Talu and Myers, 2001),
CO2CO2
=k
B
= 235:9 K
and
CO2silica
=k
B
= 147:1 K. In addition, since we also compute adsorption isotherm of nitrogen
in the same core sample, we use the LJ parameters (Ravikovitch et al., 2001; Talu and Myers,
2001),
N2N2
=k
B
= 94:45 K, and
N2silica
=k
B
= 147:3 K.
In principle, the parameter
i
can also be determined by MD simulations (Bakhshian and
Sahimi, 2017). But, here, we treat it as an adjustable parameter in order to t the model to
the data. Thus,
i
, the model's single tting parameter, was estimated by tting the calculated
sorption isotherms to our experimental data for CO
2
and N
2
adsorption in the core samples at
195 K and 77.29 K, respectively. The optimal values of
=
i
=
ff
were determined to be40
and30 for CO
2
and N
2
, respectively. Negative values of
imply that the porous medium swells
as a result of gas adsorption.
As described previously, we use the mean-eld density-functional theory (DFT) to compute
the contribution of the gas to the total energy of the system. As in any DFT formulation, one
43
must estimate the maximum adsorption capacity, n
max
v
p
, with n
max
and v
p
being, respectively,
the maximum density of the adsorbed phase, and the sorption pore volume. Their estimates are
necessary because the gas density and the adsorbed amounts predicted by the DFT are in units
of lattice site occupancy and, thus, must be converted to physical units. The minimum value of
n
max
is equal to reciprocal of the van der Waals co-volume, which for CO
2
is 1.03 g/cm
3
and
agrees with the experimental data (Hocker et al., 2003). The maximum volume of n
max
is the
number density of a close packing of gas particles, assumed to be spherical, which is
p
2=
3
with
being the molecular diameter of the gas, which for CO
2
is 1.74 g/cm
3
. Since the minimum
theoretical value of n
max
agrees with the experimental data, we used it in our calculations. The
corresponding value for N
2
is 0.723 g/cm
3
. The sorption pore volumev
p
is that of the core sample,
computed based on the experimental values of the gas density and the sample's porosity whose
sorption isotherm we have measured.
3.5.2 The Permeability of the Core Sample
Measurement of the gas permeability K of the core sample yielded, K = 325 millidarcy (mD).
As described in the previous section, we computed the permeability of the core sample using its
image together with the lattice-Boltzmann method. The result turned out to be, K = 368 mD,
which diers from the experimental result by about 13 percent.
3.5.3 Adsorption and Volumetric Strain Isotherms
The absolute and excess adsorption were calculated by the equations given in Chapter 2, Section
2.4.1. The computations were carried out as a function of the pressure P . The bulk density
b
used in the DFT is in units of lattice site occupancy and is related to the physical bulk density
p
through,
b
=
p
=n
max
. To estimate
p
, one may use either the experimental data or an
equation of state. Since the Peng-Robinson equation of state (Peng and Robinson, 1976) provides
44
(Lopez-Echeverry et al., 2017) accurate predictions for the density of gases over a wide range of
pressure, we utilized it for calculating
p
. The bulk chemical potential
b
is given by the standard
mean-eld equation (Monson, 2012)
b
=k
B
T ln
b
1
b
Z
ff
b
; (3.6)
where Z is the coordination number. In the model Z = 4, as we used tetrahedral grid blocks.
A series of bulk pressuresP were selected and the corresponding bulk density
p
were calculated
using the Peng-Robinson equation of state. Then, the corresponding bulk DFT densities,
b
=
p
=n
max
were computed, and Eq. 3.6 was used to compute the bulk chemical potential
b
. The
results were then used to compute the sorption isotherms along with the mechanical response of
the porous solid. The bulk pressure does, of course, exert an equivalent external force on the
boundaries of the porous sample, which is utilized in the calculation of the volumetric strain via
the nite-element formulation. Figure 3.1 presents the 3D image of the porous medium and its
grid structure used in the computations. The shear modulus G and the Poisson's ratio of the
solid matrix of the core sample were taken to be 2.73 GPa and 0.26, respectively, which are the
typical experimental values reported for sandstone.
The computed sorption isotherms are shown in Figure 3.2, where they are also compared with
our experimental data. It is clear that the data and the computed isotherms agree closely. Also
shown are the computed desorption isotherms in the same core sample. The dierence between
adsorption and desorption isotherms for both gases is small, although it is a bit larger in the case
of N
2
. According to the classication of the International Union of Pure and Applied Chemistry
(Sing et al., 1985) (IUPAC), the sorption isotherm of N
2
is of Type III, which corresponds to
a case in which attractive adsorbent-adsorbate interactions are weak, but molecular interactions
between the adsorbates themselves are strong. On the other hand, although Figure 3.2 does
not indicate it because the calculations and measurements in Figure 3.2 are restricted to low
45
Figure 3.1: Image of the sandstone sample (left) and its discretization with adaptive tetrahedral
mesh. Green represents that pore space, while blue indicates the solid skeleton. The image
was segmented and discretized using Seg3D 2.4.0 (CIBC, 2016a) and Cleaver 2.2 (CIBC, 2016b),
respectively.
pressures, in the IUPAC classication the CO
2
sorption isotherm corresponds to Type I, which
occurs when adsorption is limited to at most a few molecular layers, so that as pressure increases,
the amount adsorbed reaches a constant value and does not change at still higher pressures. This
is demonstrated shortly.
The volumetric response of the solid matrix of the porous medium represents its adsorption-
induced swelling and elastic compression arising due to the
uid pressure. Interestingly, the
volumetric strain \isotherms" in the solid matrix of the core sample, calculated as the sum of the
principal strains, corresponds closely to the sorption isotherms of the two gases. This is shown in
Figure 3.3, where we present the pressure-dependence of strain during sorption of the two gases.
The deformation of the rock due to adsorption of the gases is small, because the temperature is
very low and the core sample is mostly quartz, which is very dicult to deform.
Since the focus of the work is sequestration of CO
2
in rock, hereafter we present the results only
for CO
2
. Thus, we simulated sorption of CO
2
in the core sample at relatively high temperatures
of 313 K and 350 K and a wide range of pressures that are relevant to the geological conditions
for CO
2
sequestration in Mt. Simon. Figure 3.4 presents the absolute and excess adsorption
46
Figure 3.2: Comparison of the experimental data with the computed adsorption isotherms. Also
shown are the computed desorption isotherms.
Figure 3.3: Computed dependence of the volumetric strains corresponding to the isotherms of
Figure 3.2.
47
Figure 3.4: CO
2
(a) absolute and (b) excess adsorption in the porous sample at 313 K and 350
K. No hydrostatic pressure was applied to the structure.
uptakes of CO
2
at the two temperatures, demonstrating that the isotherms are indeed of the
aforementioned Type I. The excess adsorption attains a maximum, followed by a decreasing trend
at higher pressures. The reason for the trend is that, at pressures below the maximum the gas
densities in the pore space and the bulk both increase with the same rate. At pressures higher
than the pressure at which the maxima are attained, however, the pore space is saturated and,
thus, the gas density and the amount adsorbed no longer change (Chen et al., 2016; Zhou et al.,
2007). The smaller magnitude of the maximum at the higher temperature and its shift to a higher
pressure are, of course, expected. At low pressures the dierences between the two isotherms are
negligible, consistent with previous experimental and theoretical studies of CO
2
adsorption in
shale and coal formations (Zhang et al., 2016; Tang and Ripepi, 2017).
The corresponding volumetric strain
V
and porosity change, induced by CO
2
adsorption in
the porous formation at 313 K and 350 K, are depicted in Figure 3.5.
V
is dependent upon the
temperature and is higher at lower temperatures, because the amount of gas adsorbed is higher
at lowT . The strain \isotherm" follows a pattern similar to sorption, namely, it increases rapidly
for pressures up to 10 MPa at 313 K and 15 MPa at 350 K, whereas at still higher pressures the
48
Figure 3.5: (a) Volumetric strain and (b) the porosity, corresponding to the isotherms of Figure
3.4, at 313 K and 350 K.
strain changes only sightly, because the pore space has been saturated by CO
2
. Similar trends
were reported for CO
2
adsorption in coalbeds (Chen et al., 2012; Pan and Connell, 2011). The
maximum
V
due to the adsorption is 0.7% and 0.68% for 313 K and 350 K, respectively, in line
with the experimental data for CO
2
adsorption on coalbeds (Tang and Ripepi, 2017). Note that
the pressure at which the strain reaches its maximum corresponds to that at which maximum
adsorption is obtained. On the other hand, the porosity increase in the porous formation induced
by CO
2
adsorption is only about 1.15% and 1.11% at 313 K and 350 K, respectively. The small
change in the porosity is, once again, re
ective of the fact that the core sample is largely quartz.
3.5.4 Eect of the Hydrostatic Stress
Since the optimal sites for CO
2
sequestration are located deep within an underground porous
formation, it is essential to consider the eect of external conning stress on the adsorption
and deformation. Thus, we computed adsorption isotherms and volumetric strains for conning
pressure of 50 and 100 MPa. In other word, in addition to the pressure of the gas in the bulk,
49
Figure 3.6: Absolute and excess adsorption, as well as the volumetric strain in the porous medium
under a hydrostatic pressure of 50 MPa and two temperatures.
50
the porous medium experiences an excess stress from the hydrostatic pressure. The results for
the hydrostatic pressure of 50 MPa are shown in Figure 3.6. The maximum values of the excess
adsorption at 313 K and 350 K are, respectively, 1.11 mmol/gr and 0.95 mmol/gr, while the
corresponding volumetric strains are 0.69% and 0.67%. These should be compared with the results
shown in Figure 3.5. The inset of Figure 3.6(c) indicates negative strain, or solid contraction, at
low pressures, followed by increasing positive strain, indicating swelling. Thus, at low pressures the
stress arising from the conning pressure dominates the eect of swelling induced by adsorption.
The eect is more prominent at higher temperatures. But, at high pressures, it is adsorption-
induced swelling that is more eective than the hydrostatic compression.
The results for a hydrostatic pressure of 100 MPa are shown in Figure 3.7. They demonstrate
patterns similar to those shown in Figure 3.6, with only slight changes in the numerical values
of the maximum excess adsorption and the volumetric strain, and with the dierence that the
contraction-dominated stage occurs over a broader pressure range and higher compressive strains.
To better understand the eect of the external pressure, we compare in Figure 3.8 the strain-
dependence of the absolute and excess amounts of sorbed CO
2
versus volumetric strain at 313 K
and 350 K and for three hydrostatic pressures.
3.5.5 Adsorption under Geological Conditions
Due to the variations of temperature and pressure with depth in geological formations, their eect
on the adsorbed amounts of gases cannot be neglected. To this end, we assumed that the pressure
and temperature are linear functions of the depth, using a pressure gradient of 15 MPa/km and
temperature gradient of 27.3
C/km (Tang et al., 2016). We then computed the depth-dependence
of the maximum adsorption capacity and volumetric strain for CO
2
. The results for sorption are
presented in Figure 3.9. In shallower depths, up to about 200 m, the dierence between the
excess and absolute adsorption is negligible. The shapes of the two curves are also similar to what
we have presented so far, namely, that the excess amount of adsorbed CO
2
rapidly increases to
51
Figure 3.7: Same as in Figure 3.6, except that the hydrostatic pressure is 100 MPa.
52
Figure 3.8: Strain-dependence of CO
2
excess and absolute adsorption in the porous medium under
three hydrostatic pressures and temperatures of (a,c) 313 K and (b,d) 350 K.
53
Figure 3.9: Dependence on the depth of the (maximum) absolute and excess adsorption of CO
2
with a hydrostatic pressure gradient of 15 MPa/km and temperature gradient of 27.3
C/km.
a maximum value at low depths, followed by a gradual decrease at higher depths. The depth-
dependence of the absolute adsorption is also similar to what we have presented so far: It increases
sharply at low depths, and then remains essentially unchanged in the samples that are at larger
depths. Both sets of results are consistent with available data for CO
2
sorption in shale and coal
(Yang et al., 2011; Tang et al., 2016). The results shown in Figure 3.9 and their comparison
with the results presented so far suggest that, although high temperatures in deep underground
formations adversely aect the maximum adsorption capacity, the eect of high pressures in the
same depth is strong enough to oset the unfavorable thermal eect.
3.6 Summary and Conclusions
We used the model that was described in Chapter 2 for sorption-induced deformation of porous
formations that is capable of predicting the sorption isotherms in the formations. The model
was used to study CO
2
adsorption in sandstones, along with the associated deformation that
it induces in the porous sample. The computations utilized the 3D image of the core sample
54
from Mt. Simon sandstone and, therefore, no assumption was made regarding the shapes and
sizes of the pore bodies and pore throats of the porous medium. The model's predictions for
low-temperature sorption isotherms of CO
2
and N
2
are in good agreement with the experimental
data. A study of CO
2
adsorption in the porous medium under a hydrostatic pressure indicated
stress-induced compression in the solid structure at low pressures, followed by swelling at high
pressures. Furthermore, the model is capable of predicting sorption isotherms and the resulting
deformation in deep subsurface porous formations, and indicates that in such formations the excess
adsorption quickly reaches a maximum and, then, gradually decreases, with absolute adsorption
and volumetric strain following the same trends with the depth.
55
Chapter 4
Adsorption and Deformation in Flexible Metal-Organic
Frameworks
4.1 Introduction
As already mentioned in Chapter 1, metal-organic frameworks (MOFs), also known as porous co-
ordination networks, constitute a class of porous crystalline materials that are formed by metallic
species covallently linked together through organic linkers (Alhamami et al., 2014; Mishra et al.,
2014). So far, a large range of MOF materials have been synthesized through various combination
of metal ions and organic elements (Ghysels et al., 2013; Yaghi et al., 2003; Frey, 2008; Perry et al.,
2009). But, due to the combinatorial possibilities, there are millions of potential (theoretically
possible) MOF materials. Thus, large-scale molecular simulations have been used (Wilmer et al.,
2011; Chung et al., 2014; Bae et al., 2010; G omez-Gualdr on et al., 2016) in order to evaluate their
various properties, such as porosity and pore surface areas, and assess their potential applications
to such important problems as hydrogen storage and CO
2
adsorption.
The type of metal and organic subunits greatly in
uences the properties of MOFs. These
materials have several important characteristics, including large surface area, tunable porosity, and
the ability to undergo reversible deformation, which explains why they are nding a remarkable
number of applications in gas storage, separation, and adsorption technologies (Li et al., 1999;
56
Chen et al., 2001; Furukawa et al., 2010; Mishra et al., 2013; Mason et al., 2011; Finsy et al., 2009;
Alaerts et al., 2008; Herm et al., 2011). In particular, some of MOFs exhibit structural transition
and deformation with large amplitudes, if exposed to such external stimuli as light, temperature,
electric eld, and pressure, or take part in gas adsorption or liquid absorption (Ghysels et al.,
2013; Beurroies et al., 2010; Liu et al., 2008). Among them is a subset of MOFs called MIL-53
(where MIL stands for Materials of Institut Lavoisier). They consist of corner metals (such as Al,
Cr, Ga, Fe) that are connected with benzenedicarboxylate ligands, giving rise to essentially one-
dimensional diamond-shaped pore channels (Mishra et al., 2014; Neimark et al., 2011; Bourrelly
et al., 2005). They represent
exible MOFs and, when exposed to gas adsorbates, undergo what is
usually referred to as \breathing" phenomenon, (Ghysels et al., 2013) which is a reversible
exing
or deformation of the framework as a function of the amount and type of the adsorbates. In
particular, MIL-53 structure undergoes a transformation from a large pore (lp) to a narrow pore
(np) state, and switches back to the lp state during the breathing phenomenon, when exposed to
such adsorbate molecules as CO
2
and CH
4
. This implies that multiple stable states can exist in
a single framework. Such unusual and useful characteristics of this class of MOFs have provided
the impetus for researchers to consider innovative applications for
exible MOFs.
Several studies have demonstrated step-wise adsorption isotherms in some MOFs, including
MIL-53 and MIL-88. Such a complex behavior is attributed to the occurrence of a sort of a
rst-order phase transition (see below), including breathing, or crystal-to-crystal transformation
in the structure of MOFs. In particular, it has been shown that MIL-88 may swell up to 200
percent in cell volume during
uid adsorption, while the amplitude of breathing may reach about
38 percent in MIL-53, when it adsorbs CO
2
(Kitagawa and Uemura, 2005; Serre et al., 2007;
Mellot-Draznieks et al., 2005). It is clear that such phenomena are due to the interaction of the
guest molecules with the structure of MOFs.
57
The aforementioned deformation and structural transition in
exible MOFs that are induced
by adsorption have been studied both experimentally and theoretically (Serre et al., 2007; Mellot-
Draznieks et al., 2005; Mota et al., 2017; Bon et al., 2015; Simon et al., 2017; Reichenbach et al.,
2011; Tanaka et al., 2015; Wu et al., 2016; Yang et al., 2013; Yue et al., 2015). Coudert et al.
(2008) studied gas adsorption and the structural transition that it induces in such materials
by formulating the problem based on the thermodynamics of the system, and calculated the
variation of the free energy of the system during the phase transition in the structure of the
material. Neimark et al. (2010) proposed a model of the breathing transition in MIL-53 based
on the variation of the stress during adsorption of a gas. Their model predicted coexistence of
the aforementioned np and lp phases at pressures close to one at which the breathing transition
occurs. Neimark et al. (2011) also studied adsorption-induced structural transition in MIL-53 (Cr).
They found that adsorption of guest molecules in
exible MIL-53 induces small reversible elastic
deformation of around 2 and 4 percent in, respectively, the lp and np phases. In addition, their
model predicted an irreversible plastic deformation of about 40 percent in the structure during
its transition from the np to lp phases. Another thermodynamic model for adsorption-induced
breathing transition in MIL-53 (Cr) was proposed by Ghysels et al. (2013). Using a parameterized
free energy model, they derived analytic expressions for the free energy of the guest molecules,
the host material and the host-guest interaction as a function of cell parameters.
In this chapter we use the statistical mechanical model of adsorption-induced deformation
that we described in Chapter 2 in order to study the same phenomena in
exible MOFs. The
model takes into account not only the interaction of the gas molecules with the material, but also
between the guest molecules themselves. The model is then used to study adsorption of CO
2
and
CH
4
in MIL-53(Al) at various temperatures, and the structural transitions that occur in them.
The rest of this chapter is as follows. The procedure for the numerical simulation of the
theoretical model is described, after which the results are presented and discussed. The chapter
is summarized in the last section.
58
Figure 4.1: Schematic of the MIL-53 framework with the cell parameters a, b and c. The view is
along the c axis. (a) Large pore, and (b) narrow pore structures.
4.2 Computer Simulation
The details of the model and numerical simulation method were described in Chapters 2 and 3.
Here, we describe only aspects that are specic to MOFs. The crystalline structures of the np
and lp phases with the cell parameters a, b and c are shown in Figure 4.1. The view is along the
c axis. For the large- and narrow-pore shapes the parameters were,
4
respectively, a=b = 1:3 and
= 14
, and a=b = 2:7 and = 46
. The shear modulus G and the Poisson's ratio were taken
to be, respectively, 7.5 GPa and 0.2 for the np phase, and 1.5 GPa and 0.2 for the lp phase, which
are the typical values reported (Neimark et al., 2011) for MIL-53 MOF.
Figure 4.2 presents the 3D images of the MOF structure that was discretized with adaptive
tetrahedral elements. The computational grid must be resolved enough to represent the accurate
morphology of the structure, and generate reasonable computational time and accurate adsorption
isotherms. A computational domain with 478,710 tetrahedra was used for simulating sorption and
deformation for the np phase, while the corresponding grid for the lp phase consisted of 438,940
tetrahedral elements. Both were selected based on some preliminary simulations that indicated
that the adsorption isotherms would not change with higher-resolution grid.
59
Figure 4.2: Discretized structures of the (a) lp, and (b) the np frameworks, using tetrahedral cells.
4.3 Results and Discussion
We have applied the model to sorption of CO
2
and CH
4
in
exible MOF-53(Al). In order to carry
out the computations, however, we must rst supply the parameters of the model. Thus, we rst
describe their estimation.
4.3.1 Parameter Estimation
The three parameters of the model, namely, the
uid-
uid interaction parameter
ff
, the
uid-
solid material interaction parameter
fs
, and the coupling strain energy parameter
i
, may be
estimated by various methods. For example, the three parameters may be estimated by molecular
dynamics (MD) simulations. Alternatively and more simply,
ff
and
fs
may be taken to be the
Lennard-Jones energy parameters for the interactions between, respectively, the gas molecules
themselves, and the gas and the surface of the material. In the spirit of the mean-eld DFT
that we utilize, we also estimate
ff
using the mean-eld theory that relates the bulk critical
temperature T
c
to
ff
, namely,
ff
=k
B
= 4T
c
=Z. Thus, with Z = 4 we obtain
ff
=k
B
= T
c
.
Therefore,
ff
=k
B
304:25 K and 190.82 K for, respectively CO
2
and CH
4
. We estimate the
other two parameters,
fs
and
i
, by tting the model to the available experimental data at
various temperatures, although, as we show below, the two parameters are not very sensitive to
60
Figure 4.3: Comparison of computed CO
2
adsorption isotherms with and experimental data.
temperature. Estimates of the two parameters are presented in Tables 4.1 and 4.2. Note that
negative and positive values of
i
correspond, respectively, to swelling and contraction of the
porous materials. We will come back to this point shortly.
4.3.2 Adsorption Isotherms
The amount adsorbed is calculated by the equation given in Chapter 2, Section 2.4.1. The
calculated adsorbed amount or
uid density is in terms of lattice units. In order to compare the
computed isotherms with the experimental data, we must convert the units to the physical ones.
This is done through the maximum adsorption capacity. Since this value for each phase of the
61
Figure 4.4: Comparison of computed CH
4
adsorption isotherms with the experimental data.
exible MOF that we study has already been reported (Boutin et al., 2010), we utilize it in our
computations to convert the computed isotherms to physical units.
To compute the adsorption isotherms, the bulk chemical potential
b
is increased incrementally
from a very large negative value to the saturation bulk chemical potential,Z
ff
=2. To convert
the chemical potential to pressure, the standard mean-eld equations of state of a lattice gas is
utilized, which are given by Eqs. 2.14 and 2.15.
Figure 4.3 compares the computed adsorption isotherms of CO
2
in MIL-53(Al) with the cor-
responding experimental data at the same temperatures. The isotherms at 254 K, 273 K and 298
K undergo a jump, somewhat similar to a rst-order (discontinuous) phase transition, which is
typical for Type IV isotherms (Sing et al., 1985; Sahimi, 2011) that mesoporous materials exhibit.
62
The sharp transition for intermediate values of pressure emanates from the aforementioned breath-
ing phenomena, and indicates a structural transition between the np and lp phases. When the
isotherm is of Type IV type, the np phase is the most energetically stable form, and to compute
the isotherms, the empty framework prior to loading it is in that phase. The computed isotherms
for CH
4
are presented in Figure 4.4, where they are compared with the experimental data (Boutin
et al., 2010). Similar to CO
2
, the isotherms at the three lowest temperatures are of Type IV type,
which is characterized by the large jump. On the other hand, the CO
2
isotherm at 343 K and
that of CH
4
at 250 K are of Type I, according to the classication of the International Union of
Pure and Applied Chemistry (Sing et al., 1985; Sahimi, 2011), which develop when adsorption is
limited to at most a few molecular layers that, naturally, occur at higher temperatures. In this
case it is the lp phase that constitutes the most stable phase. In the context of our theoretical
model, each phase, np or lp, is characterized by its own values of the model's three parameters.
We estimated the two tting parameters and their dependence on the temperature for both CO
2
and CH
4
; they are listed, respectively, in Tables 4.1 and 4.2 for both np and lp phases, where we
present
=
fs
=
ff
and
=
i
=
ff
. Several points are worth pointing out.
(i) The strain energy parameter
is the same for both CO
2
and CH
4
. This is expected
because
is a characteristic of the porous material.
(ii)
is negative for the lp phase, but positive for the np phase, indicating, as pointed out
earlier, contraction of the material in the former phase, but swelling in the latter phase. This
indicates the internal consistency of our theoretical model.
(iii) The absolute value of
for the lp phase is higher than that for the np phase, implying
higher volumetric strain in the lp framework.
(iv) For both CO
2
and CH
4
the parameter
for the np phase varies only weakly with temper-
ature. In fact,
= 2:85 0:15 represents an accurate average value, with an estimated variation
of only 5.2 percent. Similarly, for CH
4
the estimate
= 2:45 0:05 is a highly accurate estimate.
63
(v) For both CO
2
and CH
4
for the np phase is larger than that for the lp phase, hence
suggesting higher anity of the np framework for adsorption of both CO
2
and CH
4
. This is
in agreement with the available data (Alhamami et al., 2014; Mishra et al., 2014; Boutin et al.,
2010), and indicates once again the internal consistency of our formulation. Note also that
for
CO
2
in the np phase is higher than that for CH
4
in the same phase, indicating that CO
2
has a
stronger interaction with the framework.
(vi) For the lp phase the parameter
exhibits a larger range of variations. But, even in this
case, the variations of
for both gases is only by a factor of about 1.6.
(vii) The dierence between values of
for CO
2
and CH
4
for both phases is not large. As
Figures 4.3 and 4.4 indicate, the model reproduces the experimental data for the np phase up
to the pressures at which the transition in the isotherms commences. At higher pressures the
computed isotherms at low temperatures exhibit a plateau that corresponds to the maximum
amount of adsorbed CO
2
and CH
4
, or the saturation capacity of the np phase (Boutin et al.,
2010). The model also provides accurate isotherms in the lp phase for pressures larger than those
in the transition zone.
4.3.3 Phase Fractions
Next, we attempt to develop an understanding of the root cause of the jumps in the adsorption
isotherms that, as mentioned earlier, are representative of a structural transition. To do so we
combine our model with a phase-mixture model (Ghou and Maurin, 2010). Previous studies
Table 4.1: Dependence of the parameters
and
on temperature for CO
2
.
T (K)
(np)
(np)
(lp)
(lp)
254 3 20 2.2 -30
273 2.9 20 1.9 -30
298 2.7 20 1.6 -30
343 1.35 -30
64
Figure 4.5: Pressure-dependence of the fractions of the np and lp phases in the transition zone
during CO
2
adsorption.
(Ghysels et al., 2013; Neimark et al., 2011) indicated that there is a mixture of the np and lp
phases in the transition zone. Thus, suppose that X
np
and X
lp
are the fractions of two phases
in the pressure range in the transition zone. In other words, the idea is that in the transition
zone the system is neither pure np phase, nor lp phase, but a mixture of the two. Then, the total
amount of adsorbed CO
2
,
t
, is computed by
t
=X
np
model
np
+X
lp
model
lp
; (4.1)
where
model
np
and
model
lp
are the computed adsorbed amounts for the pure np and lp phases in the
same pressure ranges. Since X
np
+X
lp
= 1, Eq. 4.1 is rewritten as
t
=X
np
model
np
+
model
lp
(1X
np
): (4.2)
There are two unknowns, but only one equation, Eq. 4.2. Thus, to determineX
np
, the adsorption
amount in the pure np and lp phases, i.e.,
model
np
and
model
lp
, are computed in the pressures range
that corresponds to the transition zone. Then, X
np
is selected in a way that
t
matches the
experimental data.
65
Figure 4.6: Pressure-dependence of the fraction of the np and lp phases in the transition zone
during adsorption of CH
4
.
Figures 4.5 and 4.6 present, respectively, the computed X
np
and X
lp
as functions of the
pressure in the transition zone for CO
2
and CH
4
, respectively. They indicates that the fractions
of the np and lp phases in the transition region in which the breathing phenomena occurs are
strong functions of the pressure. Thus, the transition is indeed a manifestation of a structural
transition in the MOF.
4.3.4 Strain and Deformation
As described earlier, numerical simulation of the model yields both the sorption isotherms and
the displacement eld in the material, resulting from the elastic deformation. Figures 4.7 and 4.8
present, respectively, the volumetric strain \isotherms" for the np and lp phases upon adsorption
of CO
2
and CH
4
. The strain in the np phase at the three lowest temperatures is negative during
adsorption of both gases, which is indicative of contraction of the framework. At the highest
Table 4.2: Dependence of the parameters
and
on temperature for CH
4
.
T (K)
(np)
(np)
(lp)
(lp)
183 2.5 20 2.2 -30
196 2.4 20 2.0 -30
213 2.5 20 1.9 -30
250 1.8 -30
66
Figure 4.7: The computed volumetric strain in the MIL-53(Al) during CO
2
adsorption in the lp
and np frameworks.
temperature, 343 K for CO
2
and 250 K for CH
4
, strain is positive, which corresponds to the
expansion of the framework in the lp framework. These results are not only consistent with those
shown in Figures 4.3 and 4.4, but also with available experimental data (Boutin et al., 2010).
4.4 Summary and Conclusions
We used the model for sorption-induced deformation of porous materials described in Chapter
2 to study the same phenomena in MOF materials. The model produces all important aspects
of sorption of CO
2
and CH
4
in MIL-53(Al), including distinct features of the well-known np
67
Figure 4.8: The computed volumetric strain in the MIL-53(Al) during CH
4
adsorption in the lp
and np frameworks.
68
and lp phases, coexistence of the two phases at intermediate pressures that is related to a struc-
tural transition, and contraction and expansion of the material at low and high pressures and
temperatures.
Let us emphasize that although the model that has been developed in this thesis was used
to simulate and compute adsorption and deformation of gases in three distinct classes of porous
materials, it is completely general and may be applied to other types of materials, including those
in which deformation is induced by the interactions between the host material and guest molecules.
The important point to recall is that, although we assumed that the materials follow linearly
elastic deformation, other types of deformation may also be considered withing the framework of
the model, so long as the constitutive equation that describes the deformation is available.
69
Chapter 5
Computer Simulation of Eect of Deformation on the
Morphology and Flow Properties of Granular Porous Media
5.1 Introduction
Over the past two decades, several theoretical studies of deformable porous media have been
undertaken. Such studies may be divided into three major classes. In one class are macroscopic
theories that are based on volume (Whitaker, 1999) or ensemble averaging (Koch and Brady,
1985) of the continuum equations of transport. One writes down the conservation laws at the
microscale, supplements them with constitutive relations that are empirical or semi-empirical,
and then averages the equations at the macroscale in order to derive the governing equations at
that scale.
Mixture theories constitute the second class of models of deformable porous media. They are
based on the macroscale conservation laws, coupled with the second law of thermodynamics and
entropy inequality. But, as pointed out by Hassanizadeh and Gray (1990, 1979b,a), since entropy
inequality is not included in averaging at the microscale to arrive at the macroscale equations,
one cannot obtain relationships that link macroscopic thermodynamic variables. In addition, such
mixture theories do not utilize any information at the microscale.
70
The third class of models of deformable porous media represents a combination of the rst
two, which we refer to as the hybrid mixture theories (HMTs) (Hassanizadeh and Gray, 1990,
1979b,a). Originally developed by Hassanizadeh and Gray (1990, 1979b), the HMTs begin with
the microscopic conservation laws, average them to derive the macroscopic equations, and then
invoke the entropy inequality in order to derive the constitutive equations, such as the generalized
Darcy's law for slow
ow of
uids through a macroscopic deformable porous medium, or Fick's
law of diusion. This approach has been developed by Weinstein et al. (2008) and references
therein, as well as Eringen (1994), Schre
er (2002), and Zhu et al. (2010).
None of the aforementioned three classes of models can take into account the eect of the
microstructure of a porous medium as it undergoes deformation. Flow and transport processes
in any porous medium, deformable or not, are controlled by its morphology that consists of
the pore-size distribution and pore connectivity (Sahimi, 2011; Sahimi and Arbabi, 1993). For
macroscopic
ow and transport processes to occur, there has to be a sample-spanning cluster of
connected pores, which exists only if the porosity of the pore space is larger than a critical
porosity, or the percolation threshold. As a porous medium is deformed, its pore space undergoes
signicant changes. The question is then, does the pore space become better connected or loses
its connectivity as the solid matrix deforms?
The question has given rise to a fourth class of models, those that can take into account the
eect of the microstructure of a porous medium on its deformation and, hence
uid
ow and
transport. Thus, there have been some eorts to develop such models and approaches. Zhu and
Wong (1999) modeled porous media by a simple-cubic network to study
uid
ow and dilatancy
in a network of microscracks, represented by the network's bonds. Boutt and McPherson (2002)
used the discrete-element method to study cracking of rock samples, but did not study
uid
ow.
Masoud and Alexeev (2010) used a random network of laments in a polymeric matrix, similar to
percolation networks with stretching and angle-changing bonds that were studied extensively in
the 1980s and 1990s (Arbabi and Sahimi, 1993; Sahimi and Arbabi, 1993) (for a comprehensive
71
review see Sahimi (2003)) to study the diusion coecient and permeance of a deformed polymer
network. Jasinski et al. (2015) studied, both experimentally and by numerical simulations,
ow
and transport properties of Bentheim sandstones, extending the earlier works of Dautriat et al.
(2009) and Thovert and Adler (2011), and similar to the work of Arns et al. (2002, 2001). They
(Jasinski et al., 2015) discretized the linear elasticity equation and solved it numerically. To
do so, the material, represented by a cubic simulation cell, was partitioned into a number of
elementary cubes that represented either the solid matrix or the pores. The geometry of the
sample was provided by an image obtained through computed microtomography. Each elementary
cube was discretized by tetrahedra. The results of the simulations were then compared with the
experimental data. Jivkov et al. (2013) used a three-dimensional (3D) network with a coordination
number of up to 14 to study the permeability of samples of rock, but the eect of deformation
was not studied.
The focus of this chapter is on studying the eect of deformation on
uid
ow through a partic-
ular class of porous materials, namely, a packing of particles that represent unconsolidated porous
media. The advantage of a packing of particles is that much is known about its microstructure
(Sahimi, 2003; Torquato, 2002). Both regular and random packings are studied, and the deforma-
tion of the particles is assumed to be small, so that the linear theory of elasticity can be utilized.
We use a quasi-static model based on Hertz's contact theory (Johnson, 1985) in order to describe
the interaction between the particles, report on the various morphological properties of the pore
space, as well as its eective permeability as functions of the applied strain and other relevant
parameters, and discuss their implications. Our study represents a step toward gaining better
understanding of the eect of deformation, caused by a variety of factors, on
ow and transport
properties of porous media modeled by a packing of particles.
72
5.2 Generation of Random Packings
We study deformation and
ow in two regular packings, the simple-cubic (SC) and face-centered
cubic (FCC) packings, as well as random packings of spherical particles. To generate the random
packings, the particles were distributed randomly in a cubic simulation cell. The particles' size
followed a Gaussian distribution, and three particle-size distributions were utilized (see below).
The boundaries in all the directions were rigid. Initially, there may be overlaps between the
particles in the random packing. To eliminate them, the following rearrangement procedure was
used. An overlap rate between two particle was dened by (He et al., 1999),R = (r
i
+r
j
d
ij
)=(r
i
+r
j
), where d
ij
is the center-to-center distance between particles i and j so that for
overlapping particles, d
ij
< (r
i
+r
j
). For each particlei a search was conducted to identify those
that may have overlap withi. Then, for each overlapping particlej centered at R
j
a new position
was calculated by
R
(n)
j
= R
(o)
j
+
r
i
+r
j
d
ij
(R
(o)
i
R
(o)
j
); (5.1)
where R
i
and R
j
are, respectively, the position vectors of particles i and j along the ray that
connects their centers to the origin of the coordinates, and superscripts n ando indicate the new
and old vectors. If particle i overlaps with n
i
other particles, then Eq. 5.1 provides n
i
positions,
and the actual new position of particle i is given by
R
i
=
1
n
i
ni
X
j=1
R
(n)
j
; (5.2)
where R
(n)
i
is calculated by Eq. 5.1. The relaxation process was applied to every overlapped
particle. If a particle neither overlaps nor contacts others, it is moved to contact its nearest
neighbor. By iterating the process, the overlap rateR decreases gradually. To avoid any bias, the
sequence of the particle rearrangements is randomized after each iteration. The overlap rateR
73
eventually drops below a preset threshold, which we took it to be 10
3
. Typically, it took about
1500 iterations to remove practically all the overlaps.
In addition to a dense random packing that has a porosity 0:32, we also generated a loose
random packing with a porosity of 0.54. The generation of the regular FCC and SC packings is
trivial and needs no explanation.
5.3 Equations of Motion and Simulation of Deformation
The model that we use for the mechanics of the packings is similar to the discrete-element method,
except that we assume the quasi-static condition in which the displacement eld throughout the
packing is calculated particle by particle. Under equilibrium condition the governing equation for
a particle with N
c
contacts is given by
F
ext
Nc
X
c=1
F
c
= 0; (5.3)
where F
ext
and F
c
are, respectively, the external and contact forces. The contact force is given
by
F
c
=F
xc
i +F
yc
j +F
zc
k = (F
n
e
ix
)i + (F
n
e
iy
)j + (F
n
e
iz
)k; (5.4)
with e
ix
, e
iy
and e
iz
being the component of the unit vector e
i
, pointing from the center of a
particlei at the spatial position X
i
and pointing toward the center of a particlej with its center's
position at X
j
, e
i
= (X
i
X
j
)=d
ij
. Here,F
n
is the magnitude of the normal contact force between
two particles that, according to Hertz's theory, is given by (Johnson, 1985)
F
n
=
4
3
Y
e
p
R
3
n
; (5.5)
74
where
n
is the overlap between two particles in contact, given by
n
=r
i
+r
j
d
ij
; (5.6)
1
Y
e
=
1
2
i
Y
i
+
1
2
j
Y
j
; (5.7)
1
R
=
1
r
i
+
1
r
j
: (5.8)
Here, R is the relative curvature, Y
i
and Y
j
are the Young's moduli of the two particles and
i
and
j
their Poisson's ratios (if the two particles are made of dierent materials), and r
i
and r
j
are the radii of the particles in contact. The normal contact force is rst considered, after which
the contribution of the tangential (shear) force F
s
or the friction force is taken into account.
According to the Coulomb's law, one must have
F
s
F
n
; (5.9)
where is the friction coecient. All the particles in the packing must satisfy inequality 5.9. The
shear or tangential force is usually written as, F
s
=k
s
s
, with k
s
and
s
being, respectively, the
tangential stiness and displacement.
Computingk
s
that enables one to calculate the contribution of the shear force is not straight-
forward, because it depends, in general, on the history of deformation and sliding of one particle
on another. On the other hand, such contributions to the deviatoric stress tensor are believed to
be small (Lu et al., 2001) and, therefore, the normal forces contribute most to the tensor. Under
such conditions, k
s
is computed by (Mindlin, 1949; Mindlin and Deresiewicz, 1953)
k
s
=
2
2=3
e
[6(1
e
)RF
n
]
1=3
2
e
; (5.10)
75
where
e
and
e
are, respectively, the eective shear modulus and Poisson's ratio, given by,
1
e
=
2
i
i
+
2
j
j
; (5.11)
2
e
=
1
i
+
1
j
; (5.12)
with
i
being the shear modulus of particle i. Due to the dierent shear forces exerted on the
particles, k
s
varies from particle to particle and from one direction to another.
Equation 5.5 expresses a nonlinear relation between the normal contact forceF
n
and
n
, based
on whichk
n
is calculated. To simplify the problem, a linear relationship is assumed (Mindlin and
Deresiewicz, 1953) to approximate the nonlinear relation (5.5), namely, F
n
= k
n
n
, with k
n
selected so as to minimize the dierence between the assumed linear relation and Eq. 5.5. Thus,
we determine (Lu et al., 2001) the minimum of the quantity,
I
2
=
Z
0
4
3
Y
e
p
R
3
n
k
n
n
2
d
n
; (5.13)
where
is the average of all particles' overlap. Then,k
n
is determined by,@I
2
=@k
n
= 0, yielding
k
n
=
8
7
Y
e
p
R
f
: (5.14)
Note that the normal stiness k
n
is not only a function of a particle's properties, but also the
particles' overlaps. One does not have to linearize Hertz' theory and can instead solve the com-
plete nonlinear problem iteratively. Our preliminary simulations with the full nonlinear problem
indicated that for small displacements (which is assumed in this work) the linearization provides
good approximation to the actual nonlinear problem with the error being on the order of a few
percent.
76
We assume that the force-displacement relation follows the Hooke's law (Lu et al., 2001).
Under the external force applied to the packing, the particles are displaced and deformed. The
incremental displacement u
x
in the x direction is then given by
u
x
=
8
>
>
<
>
>
:
P
c
F
xc
k
n
+k
s
if
ks
kn+ks
j
P
c
F
xc
j<
q
(
P
c
jF
yc
j)
2
+ (
P
c
jF
zc
j)
2
;
P
c
F
xc
p
(
P
c
jF
yc
j)
2
+(
P
c
jF
zc
j)
2
k
n
otherwise:
(5.15)
The sums are over all the contact points. The sign must be selected in such a way that the
frictional and active forces are in the directions opposite of each other. In a similar manner, the
incremental displacements in the y and z directions are given by
u
y
=
8
>
>
>
<
>
>
>
:
P
c
F
yc
k
n
+k
s
if
ks
kn+ks
j
P
c
F
yc
j<
q
(
P
c
jF
xc
j)
2
+ (
P
c
jF
zc
j)
2
;
P
c
F
yc
q
(
P
c
jF
xc
j)
2
+ (
P
c
jF
zc
j)
2
k
n
otherwise;
(5.16)
and
u
z
=
8
>
>
>
<
>
>
>
:
P
c
F
zc
k
n
+k
s
if
ks
kn+ks
j
P
c
F
zc
j<
q
(
P
c
jF
xc
j)
2
+ (
P
c
jF
yc
j)
2
;
P
c
F
zc
q
(
P
c
jF
yc
j)
2
+ (
P
c
jF
xc
j)
2
k
n
otherwise:
(5.17)
It is clear that the displacement of the neighboring contacting particles generates a new nonequi-
librium state for a given particle. Thus, the calculations utilized an iterative process in which the
forces exerted on a given particle were gradually released until the net force on each particle in
the packing was smaller than a threshold. Typically, 150-250 iterations were needed to reach the
equilibrium state.
Oedometric compression tests were simulated with both dense and loose packings by imposing
various modes of vertical compressive strain on the top boundary. During each compression cycle
the bottom face and the lateral boundaries of the system were held xed. For each strain mode,
77
50 100 150 200 250 300
0
10
20
30
40
50
60
Grid Size N
Error (%)
50 100 150 200 250 300
0
5
10
15
20
25
30
Grid Size N
Error (%)
Dense
Non−Dense
Figure 5.1: The error,jK
N
K
250
j=K
250
in calculating the eective permeability K
N
of a grid of
size NNN relative to K for N = 250.
simulations were carried out to determine the nal conguration of the packing. The applied
strain was increased up to 3 and 30 percent for the dense and nondense packings, respectively.
Once the particles reached their equilibrium state under compression,
uid
ow was simulated in
the deformed pore space in order to compute the eective permeability and other properties of
the pore space as functions of the applied strain. The eective permeability was calculated by
the lattice-Boltzmann method, which we describe brie
y in the next section.
5.4 Calculation of the Eective Permeability
The eective permeabilities of the three types of packings were calculated using a single-relaxation
time lattice-Boltzmann (LB) simulator to simulate
uid
ow in the pore space, and invoking
Darcy's law to compute the permeability. The
ow was assumed to be in the x direction (see Fig.
5.2). The no-slip boundary condition was imposed on the
uid-solid interface, as well as the rigid
outer walls of the system (in the y and z directions) by using the standard bounce-back method,
and a constant pressure gradient was applied in the
ow direction. The computational domain
for the
ow calculations consisted ofNNN grid points, which must be resolved enough that
the eective permeability would not change if a more resolved grid is used. Thus, we carried out
a series of preliminary simulations with N = 100; 150; 200; 225 and 250 in order to determine
78
Figure 5.2: The three particle-size distribution (PSD) used in the for random packings, along
with an example of a random packing before and after compression.
an adequate resolution for the computational grid. The results are shown in Fig. 5.1, where the
relative error is dened by, error =jK
N
K
250
j=K
250
, with K
N
being the eective permeability
computed with a grid of size N
3
, and K
250
representing the calculated permeability with the
most resolved grid of linear sizeN = 250. Figure 5.1 indicates that the relative error between the
results for N = 200 and N = 250 is essentially zero. Thus, all the results presented below were
obtained with a grid of size 200
3
.
The LB simulator uses an iterative sequence of propagations and collisions of ctitious particles
at the grid points of the discretized domain. The mass density of each grid point that does not
belong to the solid matrix is dened by a set of scalar particle distribution functions f(x;t), each
of which is related to a lattice velocity unit vector. At each time step, the particle distribution
functions are shifted to neighboring grid points according to their unit velocity vector through
the propagation step, also called streaming. Following each streaming step, a collision operation
is carried out to update each distribution. The simplest form of the collision operator is based
on the Bhatnagar-Gross-Krook (BGK) approximation (Bhatnagar et al., 1954). The streaming
79
and collision of the ctitious particles are described by the following equation for the particle
distribution function f
i
(Guo et al., 2000):
f
i
(x + e
i
4t;t +4t)f
i
(x;t) =
1
R
[f
i
(x;t)f
eq
i
(x;t)]; (5.18)
where t represents the time step in lattice units,
R
is the dimensionless relaxation time, and e
i
are the velocity basis vectors. The left side of Eq. 5.18 represents the streaming step by which the
particle distribution functions are shifted based on the velocities, while the right side describes
the collision operation. The equilibrium distribution f
eq
i
is given by (Qian et al., 1992),
f
eq
i
=!
i
1 +
3
c
2
(e
i
v) +
9
2c
4
(e
i
v)
2
3
2
v
2
; (5.19)
where c = x=t is the lattice speed, x is the grid spacing, v is the velocity eld, the
uid's
density, and!
i
are the weight coecients. We use the standard D3Q19 LB model for 3D modeling
(Qian et al., 1992), for which !
i
= 1=18 for i = 1; 2;; 6, !
i
= 1=36 for i = 7; 8;; 18 and
!
i
= 1=3 fori = 19. The
uid density is given by (Chen and Doolen, 2003), =
P
i
f
i
, the velocity
by, v =
P
i
e
i
f
i
, and pressure by, P =c
2
s
, where c
s
denotes the speed of sound in lattice units
(equal to c=
p
3). The kinematic viscosity is the given by, = (
R
0:5)c
2
s
. Low-Mach number
Ma
ow is simulated with Ma =v=c
s
1. We used a relaxation parameter
R
= 1. Stewart et al.
(2006) showed that if the relaxation parameter deviates from unity, the resulting permeabilities
are not accurate.
Thus, after the porous medium attained its equilibrium state for any given strain, the LB
model of
uid
ow was used in combination with the Darcy's law
hvi =
K
rP =
K
@P
@x
; (5.20)
80
0 0.5 1 1.5 2 2.5
0
0.05
0.1
0.15
0.2
0.25
Pore Size (mm)
Distribution
Before Deformation
After Deformation
0 0.5 1 1.5 2 2.5
0
0.05
0.1
0.15
0.2
0.25
Pore Size (mm)
Distribution
Before Deformation
After Deformation
0 0.5 1 1.5 2 2.5 3
0
0.05
0.1
0.15
0.2
0.25
Pore Size (mm)
Distribution
Before Deformation
After Deformation
PSD2−Hard PSD3−Hard PSD1−Hard
Figure 5.3: The computed pore-size distributions of the dense packings of hard particles with the
three particle-size distributions shown in Fig. 5.2.
in order to calculate the eective permeability K of the deformed porous medium for the corre-
sponding strain, wherehvi denotes the volume-averaged
ow velocity, and is the
uid's dynamic
viscosity. As the Darcy's law is accurate only for creeping
ows, all the simulations were carried
out for Reynolds numbers Re< 1.
We carried out extensive simulations for packings made of two types of particles, soft and hard.
The Young's modulus, shear modulus and Poisson's ratio of the soft particles were 2.7 GPa, 0.97
GPa and 0.42, whereas those for the hard particles were 110 GPa, 44 GPa and 0.24. These are in
the range of properties reported for clay particles. Typically, we simulated a packing with about
1500 - 2000 particles. The cubic domain of the packings had dimensions of 8.5 cm, 7.5 cm, 8.14
cm and 8 cm for nondense random packing, dense random packing, and SC and FCC packings,
respectively. The sizes of the particles in the random packing were not equal, but followed a
Gaussian distribution. Three particle-size distributions were utilized with an average of 3.7 mm
and standard deviations of 0.2 mm, 0.3 mm and 0.5 mm. The size of the particles in the regular
packings was 3.7 mm. Figure 5.2 presents the three particle-size distributions that we utilized in
the simulations of the random packings, along with a typical conguration of a deformed packing.
In the gures described and discussed below, we refer to the three particle-size distributions as
PSD1, PSD2, and PSD3. Unless specied otherwise, in all cases the friction factor was assumed
to be 0.03.
81
0 2 4 6 8
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
0 2 4 6 8
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
0 2 4 6 8 10
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
PSD2−Hard PSD3−Hard
PSD1−Hard
Figure 5.4: The computed pore-length distributions of the dense packings of hard particles with
the three particle-size distributions shown in Fig. 5.2.
5.5 Results and Discussion
In a dense regular packing of particles with equal sizes, the particles move only in the loading
direction under the oedometric compression. As a result, the strain eld in the packing is almost
uniform everywhere. The entire porous medium deforms as a continuous homogeneous material.
In addition, the identical spheres rearrange themselves in a crystallike formation when an external
force is applied to the packing. In other words, the particles move homogeneously due to the
application of the external force. Consequently, the deformation of the regular packing of equal-
size particles is due mainly to the local deformation at the particle contact points, rather than
being due to particles' relocation.
In the case of random packings, however, the particles move in various directions when an
external force is applied. The strain eld is not uniform even if the external force is, as in the
oedometric compression that we simulate. This eect is particularly signicant for particles that
are located close to the system's boundary. Thus, the deformation of the particles in the random
packing is due to their rearrangement, as well as the local deformation at the contact points. In
what follows we present the results and discuss their implications.
82
5.5.1 Distributions of the Pores' Size and Length
One of the most interesting properties of the packings of particles as they undergo deformation
is the change in their pore-size distribution (PSD) and pore-length distribution (PLD), which
directly aect their eective permeability and porosity. Thus, we rst calculated the PSD and
PLD of the pore space of the random packings. To do so, a well-known result (Kerstein, 1983) was
used to map the pore space between the particles onto an equivalent 3D Voronoi network in which
the bonds of the network, representing the
ow channels between the particles, are the edges of
the Voronoi polyhedra. To carry out the mapping, we rst used (Bryant et al., 1993; Dadvar and
Sahimi, 2003) the Delaunay tessellation to divide the packing into cells of four nearest neighbor
particles, with each cell being a tetrahedron with its four vertices located at the particles' centers
in the cell. Access to this void region is through any of the four surfaces that are created by three
particles in the cell. Since the vertices of the tetrahedron are the same as the particles' centers,
the void projected onto each of the four faces indicates the smallest cross section perpendicular
to the
ow eld that would be accessible for
uid
ow into the central void space of the cell. Due
to the random structure of the packing, the sizes and shapes of such cross sections vary in space.
Delaunay tessellation is the geometric dual of the Voronoi tessellation (Watson, 1981). The
Voronoi cells are polyhedra, each of which contains exactly one particle of the packing, such that
any point within every polyhedron is closer to that particle's center than the center of any other
particle in the packing. Hence, the Delaunay tessellation is used to construct a network that
represents the solid phase of the random packing, whereas the Voronoi tessellation is utilized for
constructing the pore network equivalent of the void space. The vertices of the Voronoi polyhedra
are connected by the edges of the polyhedra that represent the pore space.
To construct the Voronoi network, the positions of network's vertices are identied and, given
the four particles that surround each vertex, the way the vertices are connected to each other are
determined. For each particle, referred to as the main particle, all the possible nearest-neighbor
83
0 0.5 1 1.5 2 2.5
0
0.05
0.1
0.15
0.2
0.25
Pore Size (mm)
Distribution
Before Deformation
After Deformation
0 0.5 1 1.5 2 2.5
0
0.05
0.1
0.15
0.2
0.25
Pore Size (mm)
Distribution
Before Deformation
After Deformation
0 0.5 1 1.5 2 2.5 3
0
0.05
0.1
0.15
0.2
0.25
Pore Size (mm)
Distribution
Before Deformation
After Deformation
PSD3−Soft PSD1−Soft PSD2−Soft
Figure 5.5: The computed pore-size distributions for the dense packings of soft particles and the
three particle-size distributions shown in Fig. 5.2.
groups of particles are considered. Then, for each group a point is located that has the same
distance from the four particles' centers, which will be a vertex of the Voronoi network if and only
if its distance from the center of the main particle is equal to, or less than, its distance from the
other particles' centers in the packing. Otherwise, such a group of four particles cannot generate
a Delaunay cell, and that point cannot be a vertex in the Voronoi network. The procedure is
repeated for all the particles in the packing, so that the spatial coordinates of all the Voronoi
vertices around each main particle, as well as the groups of four particles that generate each
Delaunay cell, are determined.
After constructing the Voronoi network, we determined the eective length and radius of each
bond (pore throat) in the network. As pointed out earlier, the cross sections of the void space
representing the
ow channels vary along their length. The length of each bond of the Voronoi
network was calculated by using the spatial coordination of its two end sites. To calculate the
eective radius of each pore throat, we dene a hydraulic diameter D
h
for each bond by the
classical relation,
D
h
= 2
ow cross section
wetted perimeter
:
Simple geometrical analysis, together with the spatial coordinates of the particles' centers in the
random yield the
ow cross sections and the wetted perimeters.
84
0 2 4 6 8
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
0 2 4 6 8
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Pore Length (mm)
Distribution
Before Deformation
After Deformation
0 2 4 6 8 10
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
PSD1−Soft PSD2−Soft PSD3−Soft
Figure 5.6: The computed pore-length distributions for the dense packings of soft particles and
the three particle-size distributions shown in Fig. 5.2.
In what follows when we refer to the PSDs and the PLDs after deformation, we mean those at
the maximum applied strain. Figure 5.3 presents the calculated PSDs of the three dense packings
of hard particles for which the particle-size distributions are shown in Fig. 5.2. The results are for
before and after deformation. Only the packing with the broadest particle-size distribution, the
PSD3, exhibits some change in its pore-size distribution after deformation. This is understandable
as the particles are hard, and the packings are dense. Figure 5.4 shows the corresponding PLDs.
In this case the dierences in the PLDs before and after deformation are relatively signicant,
since the hard particles may not deform much, but can move signicantly under compression.
Figure 5.5 shows the computed PSDs of the three dense packings of soft particles, before
and after deformation, with the particle-size distributions shown in Fig. 5.2. Once again, only
the packing with the broadest particle-size distribution exhibits some change in its PSD after
deformation. This indicates that it is the high density of the particles, rather than the rigidity
of the particles, that is responsible for the changes in the PSDs. The corresponding PLDs are
presented in Fig. 5.6. The magnitude of the dierences in the PLDs before and after deformation
is comparable to those for the packings of the hard particles, hence conrming our assertion that
it is the particle density that plays the most important role in the deformation process that we
have simulated.
85
0 0.5 1 1.5 2 2.5 3
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Pore Size (mm)
Distribution
Before Deformation
After Deformation
0 0.5 1 1.5 2 2.5 3 3.5
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
Pore Size (mm)
Distribution
Before Deformation
After Deformation
0 1 2 3 4
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Pore Size (mm)
Distribution
Before Deformation
After Deformation
PSD3−Hard PSD1−Hard
PSD2−Hard
Figure 5.7: The computed pore-size distributions of random packings of hard particles with an
initial porosity 0.54 and the three particle-size distributions shown in Fig. 5.2.
0 5 10 15
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
0 5 10 15
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
0 5 10 15
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
Before Deformation
After Deformation
Before Deformation
After Deformation
PSD2−Hard PSD1−Hard
PSD3−Hard
Figure 5.8: The computed pore-length distributions of the random packings of hard particles with
an initial porosity 0.54 and the three particle-size distributions shown in Fig. 5.2.
86
0 0.5 1 1.5 2 2.5 3
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
Pore Size (mm)
Distribution
0 0.5 1 1.5 2 2.5 3 3.5
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Pore Size (mm)
Distribution
0 1 2 3 4
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
Pore Size (mm)
Distribution
Before Deformation
Before Deformation
Before Deformation
Before Deformation
Before Deformation
Before Deformation
PSD2−Soft PSD1−Soft PSD3−Soft
Figure 5.9: Same as in Fig. 5.7, but for random packings of soft particles.
0 5 10 15
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
0 5 10 15
0
0.01
0.02
0.03
0.04
0.05
0.06
Pore Length (mm)
Distribution
0 5 10 15
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Pore Length (mm)
Distribution
Before Deformation
After Deformation
Before Deformation
After Deformation
Before Deformation
After Deformation
PSD3−Soft PSD2−Soft
PSD1−Soft
Figure 5.10: Same as in Fig. 5.8, but for random packings of soft particles.
Next, we present the results for a random packing with a porosity of 0.54. Figure 5.7 depicts
the computed PSDs of the three random packings of hard particles, before and after deformation.
Relative to the dense packings, the changes in the PSDs are more pronounced. The changes in
the PLDs are even clearer; see Fig. 5.8. This is understandable because the higher porosity of
the packings allows the particles to move farther after the external strain is applied.
Similarly, if the random packings with a porosity of 0.54 consist of soft particles, then, the
changes in their PSDs and PLDs are even larger after deformation, as shown, respectively, in Figs.
5.9 and 5.10. Note that in all the cases, the deformation gives rise to a larger number of shorter
pore throats than those in the packings before the deformation. The signicant changes in the
distribution of the pore lengths also aect the permeability of the packings after deformation.
This will be discussed shortly.
87
Table 5.1: Numerical values of the exponent m and the prefactor n [see Eq. 5.21] for various
types of packings, particle-size distributions, and type of particles (soft and hard).
Loose random and regular
Packing m n
PSD1-hard 2.00 2.24
PSD2-hard 1.89 2.00
PSD3-hard 1.83 1.90
PSD1-soft 1.57 1.45
PSD2-soft 1.47 1.40
PSD3-soft 1.43 1.39
SC-hard 1.28 1.13
SC-soft 1.18 1.12
FCC-hard 1.54 5.92
FCC-soft 1.35 4.96
Dense
PSD1-hard 1.0 1.56
PSD2-hard 1.0 1.60
PSD3-hard 1.0 1.61
PSD1-soft 1.0 1.89
PSD2-soft 1.0 2.10
PSD3-soft 1.0 2.18
5.5.2 Porosity
The eect of deformation on the porosity is also important. The change in the porosity does
depend, of course, on the stiness of the particles and their size distribution. Figure 5.11 presents
the
uctuations in the porosities of the dense random packings with the particle-size distribution
PSD3 in two directions, namely, in the
ow direction x, and in the z direction along which the
external strain was applied. They represent the porosities in the planes perpendicular to the given
directions, and indicate that the deformation propagates well throughout the packings. Similar
results were obtained for other types of packings.
Figure 5.12 presents the strain-dependence of the porosity of the dense packings with the two
types of particles, soft and hard, and their three size distributions shown in Fig. 5.2. In this gure
and in the following the porosity is normalized with respect to its value before deformation. In
the case of a narrow particle-size distribution the porosity is reduced by only about 5 percent at
88
0 50 100 150 200 250
−0.035
−0.03
−0.025
−0.02
−0.015
−0.01
−0.005
z (lattice units)
Δφ (z)
0 50 100 150 200 250
−0.05
−0.04
−0.03
−0.02
−0.01
0
x (lattice units)
Δφ (x)
Figure 5.11: Fluctuations in the porosity in the strain (z) and
ow (x) directions.
0 0.005 0.01 0.015 0.02 0.025 0.03
0.93
0.94
0.95
0.96
0.97
0.98
0.99
1
Strain
Normalized Porosity
Soft
Hard
0 0.005 0.01 0.015 0.02 0.025 0.03
0.93
0.94
0.95
0.96
0.97
0.98
0.99
1
Strain
Normalized Porosity
Soft
Hard
0 0.005 0.01 0.015 0.02 0.025 0.03
0.93
0.94
0.95
0.96
0.97
0.98
0.99
1
Strain
Normalized Porosity
Soft
Hard
PSD3 PSD2 PSD1
Figure 5.12: Dependence of porosity of the dense packings on the external strain. The particle-size
distributions are shown in Fig. 5.2, and the porosities have been normalized by their initial values
before deformation.
the maximum strain applied, with the change being even smaller if the particles are hard. The
reduction in the porosity rises to about 7 percent with the widest particle-size distribution (PSD3
in the gure). As expected, the change in the porosity is larger when the particles are soft. Similar
behavior is obtained if we use regular packings of particles.
The reduction in the porosity is, however, signicantly higher in a random packing with an
initial porosity of 0.54. Figure 5.13 presents the results for the two types of particles and the
three particle-size distribution of Fig. 5.2. Even in the packing of the hard particles the porosity
at the maximum strain has been reduced by as much as 20 percent. The reduction in the porosity
89
0 0.05 0.1 0.15 0.2 0.25 0.3
0.7
0.75
0.8
0.85
0.9
0.95
1
Strain
Normalized Porosity
Soft
Hard
0 0.05 0.1 0.15 0.2 0.25 0.3
0.7
0.75
0.8
0.85
0.9
0.95
1
Strain
Normalized Porosity
Soft
Hard
0 0.05 0.1 0.15 0.2 0.25 0.3
0.7
0.75
0.8
0.85
0.9
0.95
1
Strain
Normalized Porosity
Soft
Hard
0.01 0.1 1
−1
−0.1
−0.01
−0.001
Strain
Φ/Φ
0
−1
0.01 0.1 1
−1
−0.1
−0.01
−0.001
Strain
Φ/Φ
0
−1
0.01 0.1 1
−1
−0.1
−0.01
−0.001
Strain
Φ/Φ
0
−1
PSD1 PSD2 PSD3
Figure 5.13: Same as in Fig. 5.12, but for loose packing with an initial porosity of 0.54. The
insets show the logarithmic plots of the same.
0 0.05 0.1 0.15 0.2 0.25 0.3
0.7
0.75
0.8
0.85
0.9
0.95
1
Strain
Normalized Porosity
Soft
Hard
0 0.05 0.1 0.15 0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Strain
Normalized Porosity
Soft
Hard
0.01 0.1 1
−1
−0.1
−0.01
Strain
Φ/Φ
0
−1
0.01 0.1 1
−1
−0.1
−0.01
Strain
Φ/Φ
0
−1
SC FCC
Figure 5.14: Same as in Fig. 5.12, but for the regular packings. The insets show the logarithmic
plots of the same.
is even larger in the case of regular packings; see Fig. 5.14 where we show the results for the SC
and FCC lattices. This is perhaps expected as the regular arrangement of the particles is more
amenable to deformation and rearrangement than those of a random packing of the particles.
We nd that in the case of nondense packing, as well as the regular ones, the changes in the
porosity follow a power law in the applied strain S,
0
1 =nS
m
; (5.21)
90
0 0.005 0.01 0.015 0.02 0.025 0.03
6
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
Strain
C
0 0.05 0.1 0.15 0.2 0.25 0.3
2
2.5
3
3.5
4
4.5
5
5.5
Strain
C
PSD3
PSD2
PSD1
PSD3
PSD2
PSD1
Non−Dense Packing Dense Packing
Figure 5.15: Dependence of the average number of contacts C on the external strain in random
packings of soft particles and the particle-size distributions shown in Fig. 5.2.
where
0
is the initial porosity before deformation. Estimates of n and m are listed in Table 5.1.
Figures 5.13 and 5.14 also show the ts. The coecient n is, of course, nonuniversal. For loose
random packings of hard particles, m 1:9 0:1, whereas with soft particles, m 1:5 0:07, for
the SC packings, m 1:23 0:05, and m 1:54 and 1.35 for the FCC packings with hard and
soft particles, respectively. Thus, the exponent m is also nonuniversal among the various types
of packings. In the case of the dense packings the porosity decreases with the strain essentially
linearly, so that, m = 1. Estimates of the coecient n for the various dense packings are also
listed in Table 5.1, while Figs. 5.12 also show the ts.
5.5.3 Number of Contacts between Particles
The mean number of contacts C, dened as the number of spheres that any particle touches, is
an important characteristic of a deforming packing. Figure 5.15 presents the results for the dense
packing of soft particles, as well as the packing with a porosity of 0.54. As the external strain
increases, so also does the number of contacts. For identical strains the increases are smaller in
the case of the nondense random packing due to the larger void space in the system.
91
−2 0 2 4 6 8
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
−2 0 2 4 6 8
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
−2 0 2 4 6 8
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
PSD3 PSD2 PSD1
Figure 5.16: The distributions of the computed normalized stresses, exerted by the
uid in the
pore space of the dense random packings of soft particles. The curves represent t of the numerical
data to the log-normal distribution, Eq. 5.23.
During the initial stages of compressing the loose packings, the particles are mostly displaced
in the large void space of the packings. So, the number of the contacts between the particles
does not change rapidly, giving rise to the concave dependence of the contact on the strain seen
in Fig. 5.15. The slope of the plot is also higher at higher strain modes. In other word, the
rate of increase of the contact number is higher at larger strains. The dense packing produces
the convex dependence of the dependence of the contact number on the strain. At the initial
stages of compression, the contact number increase rapidly because the packing is dense and the
particles come into contact easily. When the strain takes on higher values, the packing becomes
very compact, leaving not enough void space for the particles to move freely and, thus, the rate
of the increase of the contact number is lower at higher strain modes.
5.5.4 Fluid Stress Distribution and Deformation
An important characteristic of
uid
ow through a deforming porous medium is the distribution
of the stresses that are exerted by the
uid inside the pore space and its evolution with the
deformation. As long as the morphology of the porous medium, the
ow regime, and the type of
the
uid (Newtonian versus non-Newtonian) do not change, the shape of the distribution would
remain the same. But, because deformation changes the morphology, any change in the stress
92
−2 0 2 4 6
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
−2 0 2 4 6 8
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
−2 0 2 4 6 8
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
PSD2 PSD1 PSD3
Figure 5.17: The distribution of the computed normalized stresses, exerted by the
uid in the
pore space of a loose packing of soft particles with an initial porosity of 0.54. The curve show the
t of the numerical data to Eq. 5.23.
distribution would also be a re
ection of the deformation. To calculate the stress distribution, we
computed the second-rank stress tensor induced by
uid
ow in the pore space, before and after
deformation. The stress tensor is given by
=
1
2
rv + (rv)
T
; (5.22)
where the notation is the same as before, and superscript T denotes the transpose operation. The
eigenvalues of the tensor were then computed, with the largest one being the most important
ow-induced stress (Porter et al., 2005; Pham et al., 2014). We then calculated the distribution
of the largest eigenvalues.
Since we wish to construct a probability distribution function (PDF) for the calculated stresses
that is hopefully applicable to other types of porous media, we normalize the stresses,
=
(hi)=, wherehi is the mean stress, and is the standard deviation. Figure 5.16 presents
the stress distribution in the dense packings of soft particles before deformation and after the
maximum strain was applied. The particles' sizes were distributed according to the three distri-
butions shown in Fig. 5.2. The shape of the distribution suggests that it may be t to a log-normal
distribution. Suppose that the mean and standard deviation of the PDF of the un-normalized
93
−2 0 2 4 6
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
−4 −2 0 2 4 6 8
0
0.1
0.2
0.3
0.4
0.5
0.6
τ
*
Distribution
Before Deformation
After Deformation
SC
FCC
Figure 5.18: The distribution of the computed normalized stresses, exerted by the
uid in the
pore space of regular packings of soft equal-size particles. (a) Simple-cubic (SC) packing at an
external strain of 0.2, and (b) face-centered cubic (FCC) packing at an external strain of 0.2.
ln(stress) are, respectively, m and, so that we denote the distribution by lnN (m;). Then, the
PDF of
is a three-parameter distribution lnN (
;m
;) given by (Pham et al., 2014)
h(
) =
1
p
2(
)
exp
(
1
2
ln(
)m
2
)
; (5.23)
where m
= m ln
, and
=hi=, with
= exp(m +
2
=2)
p
exp(
2
) 1. Note that
represents the minimum value of the normalized stress. As Fig. 5.16 indicates, Eq. 5.23
provides reasonably accurate representation of the numerical data. Given the dense structure of
the packings, the stress distribution changes somewhat more signicantly only in the pore space
with the widest particle-size distribution.
Figure 5.17 presents the same, but for the random packings with an initial porosity of 0.54.
Due to the larger void space in the packings and, therefore, more signicant deformation of the
particles, the changes in the PDF of the stresses are much larger than those in the dense packings.
Note also that the range of the stresses exerted on the pore space with the wider particle-size
distribution is broader in the loose packings than the dense ones. In addition, in both types of
the packings the stress distributions in the undeformed state are more sharply peaked. This is,
94
0 0.005 0.01 0.015 0.02 0.025 0.03
0.8
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Strain
Normalized Permeability
0 0.005 0.01 0.015 0.02 0.025 0.03
0.8
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Strain
Normalized Permeability
0 0.005 0.01 0.015 0.02 0.025 0.03
0.8
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Strain
Normalized Permeability
Soft
Hard
Soft
Hard
Soft
Hard
PSD3 PSD1 PSD2
Figure 5.19: Dependence of the permeability of dense random packings on the external strain.
The particle-size distributions are those shown in Fig. 5.2. The permeabilities are normalized by
their values before deformation.
of course, due to pore closure or signicant reduction in the pore sizes, on the one hand, and
opening up other pores, on the other hand, both caused by deformation.
In Fig. 5.18 we present the corresponding results for the FCC and SC packings. In these cases
the distributions are no longer as sharply peaked as those in the random packings. The reason
is that the regular structure of the FCC and SC packings gives rise to broader distribution of
the local
ow velocities and, therefore, broader distribution of the stresses. Although the t of
the numerical results to the long-normal distribution is not as accurate as those for the random
packings, Eq. 5.23 does provide a reasonable description of the stresses.
The fact that in all the cases the stress distribution can be accurately represented by a log-
normal PDF is signicant. The long tail of the distribution implies that there is a small but
signicant fraction of the pore space in which very large stresses exist. This, of course, has direct
implications for the solid phase as well, because if the stress in the
uid in the vicinity of the
solid surface is large, it directly aects the surface. On the other hand, the sharp and narrow
peaks of the distributions imply that a large fraction of the pore space experience small stresses,
implying that the internal solid surface in that part of the pore space also experiences small
stress. Another characteristic of log-normal distribution is their variance, which means that there
are wide
uctuations in the stresses.
95
0.29 0.3 0.31
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Porosity
Normalized Permeability
0.305 0.31 0.315 0.32 0.325 0.33 0.335
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Porosity
Normalized Permeability
0.285 0.29 0.295 0.3 0.305 0.31 0.315
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Porosity
Normalized Permeability
PSD1 PSD2 PSD3
Figure 5.20: Dependence of the permeability of dense random packings on the porosity. The
particle-size distributions are those shown in Fig. 5.2. The permeabilities are normalized by their
values before deformation.
5.5.5 Eective Permeability
Figure 5.19 presents the strain-dependence of the permeabilities of the dense packings with the
two types of particles, soft and hard, in which the size of the particles is distributed according to
the particle-size distributions of Fig. 5.2. The permeabilities are normalized with respect to their
values before deformation. The results indicate several important features:
(i) The reduction in the permeabilities is not a strong function of the elastic moduli of the
particles.
(ii) Even though the nal applied strain is only 0.03, the reductions in the permeabilities are
very signicant, whereas as Figs. 5.12 and 5.13 indicate, the corresponding reductions in the
porosities are much smaller.
(iii) The change in the permeabilities is also not a strong function of the particle-size distri-
bution.
These features are due to the dense structure of the packings. As a result of applying an
external strain to the top surface of the system, the deformation propagates signicantly toward
the lower part of the packings, resulting in the reduction in the permeability.
96
0.93 0.94 0.95 0.96 0.97 0.98 0.99 1
0.82
0.84
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Normalized Porosity
Normalized Permeability
PSD1
PSD2
PSD3
Figure 5.21: Collapse of the permeability-porosity data of Fig. 5.20 onto a universal curve.
0 0.05 0.1 0.15 0.2 0.25 0.3
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Strain
Normalized Permeability
Soft
Hard
0 0.05 0.1 0.15 0.2 0.25 0.3
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Strain
Normalized Permeability
Soft
Hard
0 0.05 0.1 0.15 0.2 0.25 0.3
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Strain
Normalized Permeability
Soft
Hard
PSD1 PSD2 PSD3
Figure 5.22: Same as in Fig. 5.19, but for random loose packings with an initial porosity of 0.54.
97
0.4 0.45 0.5 0.55
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Porosity
Normalized Permeability
0.4 0.45 0.5 0.55
0.4
0.6
0.8
1
Porosity
Normalized Permeability
0.4 0.45 0.5 0.55 0.6
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Porosity
Normalized Permeability
PSD2 PSD1 PSD3
Figure 5.23: Dependence of the permeability of the random loose packings on the porosity. The
initial porosity was 0.54. The particle-size distributions are those shown in Fig. 5.2. The perme-
abilities are normalized by their values before deformation.
Figure 5.20 presents the dependence of the normalized permeabilities on the porosities for the
dense packings. If we replot the normalized permeability data shown in Fig. 5.20 versus the
normalized porosities, all the three curves collapse onto each other. This is shown in Fig. 5.21.
The results for the strain-dependent permeability of the random loose packings with a porosity of
0.54 are presented in Fig. 5.22. Note that the nal applied strain to the packings is ten times larger
than that in the dense packings. In this case the dierences between the permeabilities of the
packings with hard and soft particles are larger at larger strains, because the larger void space of
the packings allows larger deformation of the soft particles. Given that the initial porosity is large,
the reduction in the permeability is even larger than the dense counterpart. Figure 5.23 depicts
the corresponding porosity-dependent permeabilities. Once again, the numerical data collapse
onto a single curve, if the normalized permeability is plotted versus the normalized porosity. This
is shown in Fig. 5.24. The data collapse implies that the relative reduction of the permeability
is a universal function of the relative reduction in the porosity, independent of the particle-size
distribution.
The results for the SC and FCC packings, shown in Fig. 5.25, indicate even larger reduction in
the permeabilities. Due to their dense structure, the changes in the permeability of both regular
packings are not sensitive to the elastic modulus of the particles.
98
0.7 0.75 0.8 0.85 0.9 0.95 1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Porosity
Normalized Permeability
PSD1
PSD2
PSD3
Figure 5.24: Collapse of the permeability-porosity data of Fig. 5.23 onto a universal curve.
The conclusions are that, (i) in all the case the reduction in the permeabilities vary essentially
linearly with the applied strain; (ii) while the change in the permeability is not a strong function of
the hardness of the particles, it depends strongly on the initial morphology of the pore space, and
(iii) the normalized permeability is a universal function of the normalized porosity, independent
of the particle-size distribution.
5.5.6 Eect of the Friction Coecient
As already pointed out, larger friction between the particles leads to smaller deformation in them.
Thus, we expect greater reduction in the permeability when the friction coecient is smaller. This
is borne out by the results shown in Fig. 5.26, which presents the results for the permeability and
porosity for two friction coecients in a loose packing.
5.6 Summary and Conclusions
Extensive numerical simulations were carried out to study the eect of deformation on the morpho-
logical and
ow properties of packed beds of spherical particles, when they are under mechanical
compaction. Both random and regular packings of particles were studied. The deformation leads
99
0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48
0
0.2
0.4
0.6
0.8
1
Porosity
Normalized Permeability
SC
0.1 0.15 0.2 0.25 0.3 0.35
0
0.2
0.4
0.6
0.8
1
Porosity
Normalized Permeability
FCC
Figure 5.25: Dependence on the external strain of the permeability of the simple-cubic (SC) and
face-centered cubic (FCC) packings of particles. The permeabilities are normalized by their values
before deformation.
0 0.05 0.1 0.15 0.2 0.25 0.3
0.7
0.75
0.8
0.85
0.9
0.95
1
Strain
Normalized Porosity
μ = 0.03
μ = 0.80
0 0.05 0.1 0.15 0.2 0.25 0.3
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Strain
Normalized Permeability
μ = 0.03
μ = 0.80
Figure 5.26: Dependence on the friction coecient of the permeability and porosity of a random
packing of soft particles with an initial porosity of 0.54. The particle-size distribution is PSD3,
shown in Fig. 5.2.
100
to the closure of some pores and opening up of others and, therefore, the resulting reduction in the
permeability varies from relatively small to large, depending on the morphology of the initial pore
space. The distribution of the pores' lengths, the porosity of the pore space, and the pore-size
distribution also evolve under an external strain. It is of course the changes in the morphology of
the pore space that lead to the strong reduction of the permeability.
The Kozeny-Carman (KC) equation (Kozeny, 1927; Carman, 1937),
K =
hDi
2
36C
KC
3
T
2
(1)
2
(5.24)
is often used to correlate the permeability of packed beds of particles with their porosity , where
hDi is the mean particle diameter,T is the volume-averaged tortuosity, and C
KC
is the Kozeny-
Carman shape factor that depends on the materials that the particles are made of, and accounts
for the variations in the permeability of porous media having the same porosity, but distinct
microstructures. C
KC
is usually taken to be 2.5 for packings of spherical particles. The tortuosity
can be dened in a variety of ways (Ghanbarian et al., 2013b,a). We nd, however, that if the
KC equation is to be used to correlate the permeability with the porosity of the deforming pore
space, C
KC
must depend on . This will, however, be purely empirical, as a result of which the
equation loses its predictive power.
101
Chapter 6
Eect of Flash Heating on the Frictional Strength
Evolution of Granular Media During Slip Weakening
6.1 Introduction
An important factor that contributes to rupture propagation and the ensuing instability, and
ultimately earthquakes, is the dynamic friction of faults during earthquake (Smeraglia et al.,
2017; Collettini et al., 2013). Previous studies have attempted to demonstrate reduction of fault
friction and, hence, its strength during earthquakes with high slip rates (Kitajima et al., 2011),
a process known as slip weakening (Wibberley and Shimamoto, 2005; Hirose and Shimamoto,
2005). In other words, the general thinking is that earthquakes occur when a fault weakens with
increasing slip rates (Sone and Shimamoto, 2009; Rice, 2006). Understanding the resistance to
coseismic slip on faults is an important problem in earthquake mechanics. The frictional resistance
aects the magnitude of the reduction in shear stress, as well as the extent of destruction that
can happen during earthquake (Di Toro et al., 2004). Low-slip resistance along faults decreases
their strength and facilitates the slip propagation (Smeraglia et al., 2017; Di Toro et al., 2004).
It has been reported that the weakening phenomena during fault slip are activated essentially
by thermal processes, such as thermal pressurizing of pore
uid and frictional heating (Collettini
et al., 2013; Di Toro et al., 2011; Spagnuolo et al., 2015; Kitajima et al., 2010; Mase and Smith,
102
1987; Lachenbruch, 1980). Flash or frictional heating is a phenomenon at microscopic scale in
which heat is generated at asperity contacts due to high shear rates, i.e. fast fault motion. The
heat generated at the contact points or surfaces would increase the average temperature of the
fault surface, hence causing melting of the entire fault surface. This phenomena may lead to
reduction of frictional shear strength of contacts. Originally,
ash heating was proposed in order
to explain the variation of dry friction as a result of high slip rates in metals (Molinari et al., 1999;
Lim et al., 1989; Ettles, 1986). It was Rice (2006) who developed a theory of
ash heating at
frictional asperity contacts in order to describe the weakening behavior and dependence of fault
friction on the slip rate in rocks (Rice, 2006, 2017). Several studies in seismology have considered
ash heating as a mechanism that controls fault friction at high-rate earthquake slips.
Several experimental studies have reported on rock and gouge friction under high-slip rates
that correspond to the those occurring during earthquakes. These studies indicate the presence
of frictional weakening for slip rates larger than 0.1 m/s, independent of rock or gouge type, and
is attributed mainly to the
ash heating phenomenon (Di Toro et al., 2011). For example, an
experimental study of quartz rocks revealed a decrease of friction coecient by a factor of 3 due
to the rock weakening at subseismic slip rates (Goldsby and Tullis, 2002). Another study (O
Hara
et al., 2006) demonstrated the weakening response of the coal gouge as a result of frictional heating
at high slip rates, about 1 m/s. These results indicate the presence of various phenomena during
slip experiments, including frictional heating,
uid pressurization, and slip weakening. Using a
theory of shear transformation zone, Elbanna and Carlson proposed (Elbanna and Carlson, 2014)
a model for
ash weakening in granular media, in order to investigate the shear behavior of the
media under various sliding rates and conning pressures. Their model predicted a logarithmic
dependence of steady-state shear strength on the imposed slip velocity. Moreover, a velocity-
weakening friction model, inspired based on the
ash heating concepts in granular media, was
developed by Lucas et al. (2014) to predict the friction coecient of landslides. It was found that
the proposed model was capable of reproducing the available data for landslides on Earth and
103
perhaps other planets. The authors also reported reduction of friction coecient as the sliding
velocity increased.
To gain a better understanding of frictional weakening during earthquakes, we have carried
out computer simulation of slip in granular media using the discrete-element method (DEM).
Various slip rates, ranging from 0.05 m/s to 1 m/s, were simulated. We assess the shear response
of the granular medium to sliding velocity as we vary its magnitude. Variation of the local friction
coecient in the medium and the dynamic evolution of frictional strength, which is induced by
slip weakening, are investigated. In the following sections, we rst describe the simulation method
and explain how to incorporate the
ash heating model in the granular medium and the DEM
simulations. We then present the results and discuss them.
6.2 Theoretical Model and Computer Simulation
Let us rst describe the DEM model brie
y, after which the numerical simulation procedure is
discussed.
6.2.1 Discrete-Element Modeling
The DEM was initially developed by Cundall and Strack (Cundall and Strack, 1979; Potyondy
and Cundall, 2004) for analyzing the behavior of granular assemblies by using a linear spring-
dashpot model. To account for the behavior of the rock, a packing of spherical particles is
used. This implies that in order to optimize the computational requirement and the complexity
in the contact detection step, spheres represent the discrete elements. Over very small time
steps, the disturbances caused by an external force can propagate from any element to only its
immediate neighboring elements. Thus, in the DEM method, at any given time step only the
pair-wise interaction between the neighboring particles is considered. In the present study a
104
three-dimensional (3D) DEM model is simulated in order to carry out shear tests in a granular
gouge and study the evolution of frictional strength of such a medium during slip weakening.
To tackle with the high computational requirement, the DEM algorithm is formulated in
parallel scheme using OpenMP (Dagum and Menon, 1998) . The interaction between the elements
is modeled via linear contact models. Based on the net contact force on each element, the motion
of the elements (particles) is simulated by numerically solving Newton's law of motion. A Verlet-
type time integration scheme (Fraige and Langston, 2004; Kruggel-Emden et al., 2008) is used.
Translational and rotational motion of each particle are dened through (Remy et al., 2009; Boac
et al., 2014):
m
i
dv
i
dt
=m
i
g +
X
j
(F
n
ij
+F
t
ij
); (6.1)
I
i
dw
i
dt
=
X
j
r
i
F
t
ij
; (6.2)
where m
i
, v
i
, r
i
, I
i
and w
i
are, respectively, the mass, linear velocity, radius, moment of inertia
and angular velocity of particle i, and g is the gravitational force. A contact search is done to
identify the overlaps between the elements. The normal overlap for each contact is calculated by
n
=r
i
+r
j
jx
i
x
j
j; (6.3)
wherer
i
andr
j
are the radius of particlesi andj, and x
i
is the position vector for particlei. The
detection of contact overlaps is followed by calculation of normal and tangential contact forces,
F
n
and F
t
, which are computed using linear-elastic contact model (Cundall and Strack, 1979):
F
n
=k
n
3=2
n
; (6.4)
F
t
=k
t
t
: (6.5)
105
Here, k
n
and k
t
are the normal and tangential stiness, respectively, and
t
is the tangential
cumulative shear displacement. k
n
is given by (Potyondy and Cundall, 2004):
k
n
= 4r
Y
(6.6)
with Y
being the equivalent Young's modulus, and r
the equivalent radius of the contacting
particles. Equivalent properties for contacting particles i and j are given by (Boac et al., 2014):
R
=
1
r
i
+
1
r
j
1
; (6.7)
Y
=
1
i
2
Y
i
+
1
j
2
Y
j
1
; (6.8)
with being the Poisson's ratio.
The tangential displacement
t
is calculated in an incremental fashion by temporal integration
of the tangential component of the relative velocity at the contact point (Remy et al., 2009; Boac
et al., 2014):
t
=
Z
v
t
dt; (6.9)
where v
t
is the tangential component of the contact velocity, dened by (Boac et al., 2014):
v
t
= (v
i
v
j
)t +!
i
r
i
+!
j
r
j
(6.10)
wheret is the tangential component of the unit vector connecting the contacting particles' centers,
and !
i
and !
j
are the angular velocity of particle i and j. The tangential force is limited by
Coulomb's law of friction (Potyondy and Cundall, 2004):
F
t
s
jF
n
j; (6.11)
106
Figure 6.1: Random packing of spheres used in the shear simulation. The upper-layer particles
move with velocity V
s
in the x direction, while the lower-level particles are held xed.
Where
s
is the static friction coecient. According to Coulomb's law, sliding occurs at the
contact points when the tangential force exceeds the static friction. The time step must be small
enough to maintain the assumption of constant translational and rotationl accelerations, and
attain stable and quasi-static simulation. The appropriate number of time steps N, which is
proportional to the computational time, is estimated by (Ergenzinger et al., 2011)
N
p
Y
r
; (6.12)
where r denotes the average radius of the particles in the packing, and is the density of the
particles. The input parameters for the DEM simulation are listed in Table 6.1.
6.2.2 Friction Weakening by Flash Heating
We assume thatT is the average temperature of the fault surface. As sliding occurs and asperity
contact forms,
ash heating causes the temperature T
a
at the asperity contact to increase to
higher values for a high-speed sliding. Thus, as the contact slides, its temperature rises during
its brief lifetime, = D
a
=V
s
, where D
a
is the contact diameter and V
s
is the slip rate. The
temperature increase is essentially due to the generated heat localized at the asperity contact at
the rate of
c
V
s
, where
c
is the frictional strength of the asperity contact. It is assumed that
c
107
is nearly constant and does not change with variations of T
a
. But, whenT
a
reaches values higher
than a weakening temperature T
w
, the frictional strength of the asperity contact decreases again
to a weakened value
w
c
. For a large slip rateV
s
, the generated heat at the asperity interface
does not have sucient time to propagate quickly, and the contact point reaches the weakened
temperature T
w
at a time
w
< . By assuming a 1D heat conduction at the asperity contact,
the time
w
is given by:
w
=
V
2
s
c(T
w
T )
c
2
; (6.13)
where is the thermal diusivity, and c is the volumetric heat capacity. Thus, the weakening
velocity is calculated by
w
=D
a
=V
w
, i.e.
V
w
=
D
a
c(T
w
T )
c
2
: (6.14)
V
w
is known as the critical sliding rate that, such that for velocities greater thanV
w
the weakening
occurs at the asperity contacts, whereas at velocities lower than V
w
there is no weakening in the
granular medium or rock.
Thus, based on the model, at low sliding rates, V
s
< V
w
, the contacts are strong enough
during their entire lifetime. But, when V
s
> V
w
, or in other words, if
w
< , the contacts
spend a fraction of their lifetime, V
w
=V
s
, with high strength, while the remaining fraction is in
the weakened state. By ignoring the statistical distribution of contact diameter D
a
and other
properties, the friction coecient of the fault exhibits the following dependence on the critical
slip rate during the weakening (Rice, 2006; Lucas et al., 2014; Beeler et al., 2008):
w
=
8
>
>
<
>
>
:
s
V
s
V
w
;
s
1 +V
s
=V
w
V
s
V
w
;
(6.15)
108
Table 6.1: Micromechanical and geometrical properties of the packing.
Parameter Value
Particle Count 3000
Particle radius 3 mm
Domain Dimensions 868685 mm
Density 2650 kg/m
3
Normal Stiness 3.24710
6
N/m
Tangential Stiness 3.24710
6
N/m
Poisson's ratio 0.26
Elastic Modulus 0.69 GPa
Friction Coecient 0.5
Normal Stress 10 MPa
Shear Velocity 0.05, 0.08, 0.1, 0.2, 0.5 and 1 m/s
Time Step Length 0.005 s
Thermal Diusivity 0.5 mm
2
/s
Volumetric Heat Capacity 2.7 MPa/
C
Weakening Temperature 1000
C
where
s
is the low-speed friction coecient, i.e., the static friction coecient, and
w
represents
the weakened friction coecient at slip rates higher than the weakening limit.
To study the eect of velocity-weakening friction on shear dynamics of granular media, we
incorporated the
ash weakening theory in the DEM shear simulations. To carry out the simula-
tions, a granular packing was generated using the method described in Chapter 5, Section 5.2. A
packing with dimension of 86 86 85 mm was constructed that consists of 3000 particles. The
schematic of the packing is shown in Figure 6.1. Next, as the preshear process, the packing was
compressed by applying a normal stress
n
to the top dynamic boundary. Afterwards, the shear
test was implemented by moving the top layer of the particles at a constant horizontal velocity in
the x direction, while the particles at the lower boundary were held xed at their locations. To
keep the condition of normally-compressed medium, the upper stress boundary condition needs
to be satised during the shearing simulations. In the shear simulations, the boundaries in the y
andz directions are considered periodic. The shear stress is calculated by summation of the force
components acting on the uppermost particles along the x axis, i.e., along the shear movement.
109
Figure 6.2: Shear stress-strain for constant and variable friction factors with shear velocity = 1
m/s.
6.3 Results
We carried out shear-test simulations at conning stress of 10 MPa and shear velocities 0.05 -
1 m/s. The shear test was carried out for two cases: (1) the local friction coecient between
the particles is constant,
s
= 0:5, and (2) the local friction coecient varies based on the
ash
heating theory as described earlier, which, in what follows, we refer to as the no-weakening and
weakening cases. The elastic modulus and Poisson's ratio of the particles were taken as 0.69 GPa
and 0.26. The normal stinessk
n
was calculated by using Eq. (6.6), and listed in Table 6.1. The
tangential stinessk
t
is taken as equal to the normal stiness. The average temperature of gouge
surface and weakening temperature are chosen as 400
C and 1000
C. The thermal diusivity
and the volumetric heat capacity (c) were selected as 0.5 mm
2
/s and 2.7 MPa/
C, respectively.
We took the frictional strength of the asperity contact,
c
, to be 1/6 of the shear modulus. Several
shear velocities were applied to the upper layer of the packing in the x direction.
Figure 6.2 represents shear stress as a function of strain associated with shear velocity of 1
m/s for weakening and no-weakening cases. As the results indicate, at small strains, there is
110
Figure 6.3: Shear stress-strain for constant and variable friction factors with shear velocity = 0.5
m/s.
Figure 6.4: Shear stress-strain for constant and variable friction factors with shear velocity = 0.2
m/s.
111
Figure 6.5: Shear stress-strain for constant and variable friction factors with shear velocity = 0.1
m/s.
a negligible dierence between the frictional stress for weakening and non-weakening cases. As
strain proceeds to higher values, however, the dierence between them becomes more signicant
until they reach a steady-state value and their dierences remain almost unchanged. The steady
shear stress for the no-weakening case is around 0.22 GPa, whereas for the weakening case is
lower, around 0.15 GPa. In other words, the steady shear strength decreases almost 32%, when
thermal weakening is taken into account.
The results associated with shear velocity of 0.5 m/s have been depicted in Figure 6.3. In
this case the steady shear stress drops from 0.11 GPa to 0.09 GPa after taking into account
ash heating, which means there is a 18% decline in the steady shear stress. As the sliding
velocity decreases to 0.2 m/s, the dierence between steady shear stress of the two cases reduces;
see Figure 6.4, while the results for smaller shear velocities are shown in Figures 6.5, 6.6 and
6.7. For low values of sliding velocities, including 0.08 and 0.05 m/s, the shear responses of the
weakening and no-weakening cases are essentially identical, and the steady-state shear strength
remains almost unchanged when
ash heating is invoked. Thus, for sliding rates higher than
0.1 m/s, the weakening eect is signicant, which in line with the previous studies that had
112
Figure 6.6: Shear stress-strain for constant and variable friction factors with shear velocity = 0.08
m/s.
Figure 6.7: Shear stress-strain for constant and variable friction factors with shear velocity = 0.05
m/s.
113
Figure 6.8: Dependence of steady shear stress on shear velocity for weakening and no-weakening.
indicated that rocks tend to weaken for at slip rates higher than 0.1 m/s (Di Toro et al., 2011).
By comparing the results, we conclude that the eect of thermal weakening in more prominent at
higher velocities. Furthermore, it is clear that as the sliding velocity decreases, the shear stress
attains the steady-state limit faster.
Figure 6.8 demonstrates how the steady shear stress depends on the shear velocity in both
weakening and no-weakening limits, indicating that as the velocity increases, the eect of thermal
weakening is more pronounced.
Histograms of weakened local frictional coecient
w
between the particles of the packing
under various sliding velocities are shown in Figure 6.9, which represent the distribution of fric-
tional coecient at dierent times during the shear test. In all the simulations, the packing was
sheared at 10 MPa normal stress. As shown in Figure 6.9(a), within the initial micro-seconds of
the process, a high percentage of the particles experience weakening as the local friction coecient
is mostly between 0.3 and 0.5, for all sliding rates. For higher sliding velocities, the drop in the
friction factor is even higher.
114
Figure 6.9: Distribution of weakened friction factor for various shear velocities and at times (all
in s)(a) 2.5; (b) 500; (c) 1000, and (d) 2150.
115
At longer times, the distribution of friction coecient has a longer tail at smaller values. This
means that over longer times the weakening state is more wide spread, penetrating more deeply
inside the packing.
6.4 Summary and Conclusions
Inspired by
ash weakening process during earthquakes, we utilized the DEM to simulate shear
test in a granular gouge. The theory of
ash heating was incorporated into the DEM simulations
in order to follow the variations of the local friction coecient between the particles in the gouge.
The shear test simulations were carried out over a range of slip rates, from 0.05 to 1 m/s. This
study is an eort to understand the velocity-dependent frictional weakening in granular media.
Variation of the local friction coecient in the medium and evolution of frictional strength, which
is induced by slip weakening, have been investigated and demonstrated.
We showed that at velocities higher than 0.1 m/s, the eect of frictional weakening on the
frictional strength - shear stress-strain behavior - of the medium is nontrivial, which is consistent
with the previous studies. The results demonstrate that frictional strength decreases due to
ash
heating, and the decrease is higher at higher sliding velocities.
116
Bibliography
Alaerts, L., Maes, M., Giebeler, L., Jacobs, P. A., Martens, J. A., Denayer, J. F. M., Kirschhock,
C. E. A., and DeVos, D. E. (2008). Selective adsorption and separation of ortho-substituted
alkylaromatics with the microporous aluminum terephthalate mil-53. J. Am. Chem. Soc.,
130:14170{14178.
Alam, A. K. M. B., Niioka, M., Fujii, Y., Fukuda, D., and Kodama, J. (2014). Eects of conning
pressure on permeability of three rock types under compression. Int. J. Rock Mech. Min. Sci.,
65:49{61.
Alhamami, M., Doan, H., and Cheng, C. H. (2014). A review on breathing behaviors of metal-
organic frameworks (mofs) for gas adsorption. Materials, 7:3198{3250.
Anggara, F., Sasaki, K., Rodrigues, S., and Sugai, Y. (2014). The eect of megascopic texture on
swelling of a low-rank coal in supercritical carbon dioxide. Int. J. Coal Geol., 125:45{56.
Arbabi, S. and Sahimi, M. (1993). Mechanics of disordered solids. i. percolation on elastic networks
with central forces. Phys. Rev. B, 47:695.
Arns, C., Knackstedt, M., Pinczewski, W., , and Garboczi, E. (2002). Computation of linear
elastic properties from microtomographic images: Methodology and agreement between theory
and experiment. Geophys., 67:1396.
Arns, C., Knackstedt, M., Pinczewski, W., and Lindquist, W. (2001). Accurate computation of
transport properties from microtomographic images. Geophys. Res. Lett., 28:3361.
Bae, Y. S., Yazaydin, A. O., and Snurr, R. Q. (2010). Evaluation of the bet method for determining
surface areas of mofs and zeolites that contain ultra-micropores. Langmuir, 26:5475{5483.
Bakhshian, S. and Sahimi, M. (2016). Computer simulation of the eect of deformation on the
morphology and
ow properties of porous media. Physical Review E, 94:042903.
Bakhshian, S. and Sahimi, M. (2017). Adsorption-induced swelling of porous media. Int. J.
Greenhouse Gas Control, 57:pp. 1{13.
Balzer, C., Cimino, R. T., Gor, G. Y., Neimark, A. V., and Reichenauer, G. (2016). Deformation
of microporous carbons during n
2
, ar, and co
2
adsorption: Insight from the density functional
theory. Langmuir, 32:8265{8274.
Balzer, C., Wildhage, T., Braxmeier, S., Reihenauer, G., and Olivier, J. P. (2011). Deformation
of porous carbons upon adsorption. Langmuir, 27:2553{2560.
Bathija, A. (2009). Elastic Properties of Clays. Ph.D. Thesis, Colorado School of Mines.
Beeler, N., Tullis, T., and Goldsby, D. (2008). Constitutive relationships and physical basis of
fault strength due to
ash heating. J. Geophys. Res., 113:B01401.
117
Benson, C., Zhai, H., and Wang, X. (1994). Estimating hydraulic conductivity of compacted clay
liners. J. Geotech. Eng., 120:368387.
Beurroies, I., Boulhout, M., Llewellyn, P. L., Kuchta, B., Frey, G., Serre, C., and Denoyel, R.
(2010). Using pressure to provoke the structural transition of metalorganic frameworks. Angew.
Chem. Int. Ed., 49:7526{7529.
Bhatnagar, P., Gross, E., and Krook, M. (1954). A model for collision processes in gases. i. small
amplitude processes in charged and neutral one-component systems. Phys. Rev., 94:511.
Billemont, P., Coasne, B., and Weireld, G. D. (2013). Adsorption of carbon dioxide, methane, and
their mixtures in porous carbons: Eect of surface chemistry, water content, and pore disorder.
Langmuir, 29:3328{3338.
Biot, M. (1941). General theory of three dimensional consolidation. J. Appl. Phys., 12:155.
Biot, M. (1956). Theory of deformation of a porous viscoelastic anisotropic solid. J. Appl. Phys.,
27:459.
Boac, J., Ambrose, R., Casada, M., Maghirang, R., and Maier, D. (2014). Applications of discrete
element method in modeling of grain postharvest operations. Food Eng. Rev., 6:128149.
Bon, V., Klein, N., Senkovska, I., Heerwig, A., Getzschmann, J., Wallacher, D., Zizak, I., Brzhezin-
skaya, M., Mueller, U., and Kaskel, S. (2015). Exceptional adsorption-induced cluster and
network deformation in
exible metal-organic framework dut-8 (ni) observed by in-situ x-ray
diraction and exafs. Phys. Chem. Chem. Phys., 17:17471{17479.
Botan, A., Rotenberg, B., Marry, V., Turq, P., and Noetinger, B. (2010). Carbon dioxide in
montmorillonite clay hydrates: thermodynamics, structure, and transport from molecular sim-
ulation. J. Phys. Chem. C, 114:14962{14969.
Bourrelly, S., Llewellyn, P. L., Serre, C., Millange, F., Loiseau, T., and Frey, G. (2005). Dif-
ferent adsorption behaviors of methane and carbon dioxide in the isotypic nanoporous metal
terephthalates mil-53and mil-47. J. Am. Chem. Soc., 127:13519{13521.
Boutin, A., Coudert, F. X., Springuel-Huet, M. A., Neimark, A. V., Ferey, G., and Fuchs, A.
(2010). The behavior of
exible mil-53(al) upon ch
4
and co
2
adsorption. J. Phys. Chem. C,
114:22237{22244.
Boutt, D. and McPherson, B. (2002). Simulation of sedimentary rock deformation: Lab-scale
model calibration and parameterization. Geophys. Res. Lett., 29:1954.
Brochard, L., Vandamme, M., and Pellenq, R. J. M. (2012). Poromechanics of microporous media.
J. Mech. Phys. Solids, 60:606 622.
Bryant, S., King, P., and Mellor, D. (1993). Network model evaluation of permeability and spatial
correlation in real random sphere packing. Transp. Porous Media, 11:53.
Busch, A., Alles, S., Gensterblum, V., Prinz, D., Dewhurst, D. N., Raven, M. D., Stanjek, H.,
and Krooss, B. M. (2008). Carbon dioxide storage potential of shales. Int. J. Greenhouse Gas
Control, 2:297{308.
Cai, B., Li, Q., Liu, G., Liu, L., Jin, T., and Shi, H. (2017). Environmental concern-based site
screening of carbon dioxide geological storage in china. Sci. Rep., 7:7598.
Carbone, A., Chiaia, B. M., Frigo, B., and T urk, C. (2010). Snow metamorphism: A fractal
approach. Phys. Rev. E, 82:036103.
118
Carman, P. (1937). ,
uid
ow through granular beds. Trans. Inst. Chem. Eng. London, 15:150.
Chareonsuppanimit, P., Mohammad, S. A., Robinson, R. L., and Gasem, K. A. M. (2012). High-
pressure adsorption of gases on shales: Measurements and modeling. Int. J. Coal Geol., 95:34.
Chen, B., Eddaoudi, M., Hyde, S. T.,
OKeee, M., and Yaghi, O. M. (2001). Interwoven metal-
organic framework on a periodic minimal surface with extra-large pores. Science, 291:1021{1023.
Chen, G., Lu, S., Zhang, J., Xue, Q., Han, T., Xue, H., Tian, S., Li, J., Xu, C., Pervukhina, M.,
and Clennell, B. (2016). Research of co
2
and n
2
adsorption behavior in k-illite slit pores by
gcmc method. Sci. Rep., 6:paper 37579.
Chen, G., Yang, J., and Liu, Z. (2012). Method for simultaneous measure of sorption and swelling
of the block coal under high gas pressure. Energy and Fuels, 26:4583 4589.
Chen, S. and Doolen, G. (2003). Lattice boltzmann method for
uid
ows. Annu. Rev. Fluid
Mech., 30:329.
Chen, T., Feng, X. T., and Pan, Z. (2015). Experimental study of swelling of organic rich shale
in methane. Int. J. Coal Geol., 150-151:64{73.
Chin, L. and Nagel, N. (2004). Modeling of subsidence and reservoir compaction under water
ood
operations. Int. J. Geomech., 1:28{34.
Cho, H., Deng, H., Miyasaka, K., Cho, Z. D. M., Neimark, A. V., Kang, J. K., Yaghi, O. M., and
Terasaki, O. (2015). Extra adsorption and adsorbate superlattice formation in metal-organic
frameworks. Nature, 527:503{507.
Chung, Y., Camp, J., Haranczyk, M., Sikora, B. J., Bury, W., Krungleviciute, V., Yildirim,
T., Farha, O. K., Sholl, D. S., and Snurr, R. Q. (2014). Computation-ready, experimental
metal.organic frameworks: A tool to enable high-throughput screening of nanoporous crystals.
Chem. Mater., 26:6185{6192.
CIBC (2016a). Seg3D: Volumetric Image Segmentation and Visualization. Scientic Computing
and Imaging Institute (SCI), Download from: http://www.seg3d.org.
CIBC (2016b). Cleaver: A MultiMaterial Tetrahedral Meshing Library and Ap-
plication. Scientic Computing and Imaging Institute (SCI), Download from:
http://www.sci.utah.edu/cibc/software.html.
Collettini, C., Viti, C., Tesei, T., and Mollo, S. (2013). Thermal decomposition along natural
carbonate faults during earthquakes. Geology, 41:927.
Coudert, F.-X., Jeroy, M., Fuchs, A. H., Boutin, A., and Mellot-Draznieks, C. (2008). Ther-
modynamics of guest-induced structural transitions in hybrid organic-inorganic frameworks. J.
Am. Chem. Soc., 130:14294{14302.
Cui, X., Bustin, R. M., and Chikatamarla, L. (2007). Adsorption-induced coal swelling and
stress, implications for methane production and acid gas sequestration into coal seams. Journal
of Geophysical Research-Solid Earth, 112:B10202.
Cundall, P. and Strack, O. (1979). A discrete numerical model for granular assemblies. Geotech-
nique, 29:4765.
Dadvar, M. and Sahimi, M. (2003). Pore network model of deactivation of immobilized glucose
isomerase in packed-bed reactors. iii: Multiscale modeling. Chem. Eng. Sci., 58:4935.
119
Dadwhal, M., Kim, T. W., Sahimi, M., and Tsotsis, T. T. (2008). The study of co
2
diusion and
adsorption on calcined layered double hydroxides: the eect of particle size. Ind. Eng. Chem.
Res., 47:6150{6157.
Dagum, L. and Menon, R. (1998). Openmp: An industry-standard api for shared-memory pro-
gramming. IEEE Comput. Sci. Eng., 5:4655.
Dang, D., Zhao, L., Lu, X., Xu, J., Sang, P., Guo, S., Zhu, H., and Guo, W. (2017). Molecular sim-
ulation of co
2
/ch
4
adsorption in brown coal: eect of oxygen-, nitrogen-, and sulfur-containing
functional groups. Appl Surf Sci, 423:33{42.
Dashtian, H., Wang, H., and Sahimi, M. (2017). Nucleation of salt crystals in clay minerals:
molecular dynamics simulation. J. Phys. Chem. Lett., 8:3166.
Dautriat, J., Gland, N., Guelard, J., Dimanov, A., and Raphanel, J. L. (2009). Axial and radial
permeability evolutions of compressed sandstones: end eects and shear-band induced perme-
ability anisotropy. Pure Appl. Geophys., 166:1037{1061.
Day, S., Fry, R., and Sakurovs, R. (2008). Swelling of australian coals in supercritical co
2
. Int. J.
Coal Geol., 74:41{52.
de Jong, S., Spiers, C., and Busch, A. (2014). Development of swelling strain in smectite clays
through exposure to carbon dioxide. Int. J. Greenhouse Gas Control, 24:149{161.
Di Toro, G., Goldsby, D., and Tullis, T. (2004). Friction falls towards zero in quartz rock as slip
velocity approaches seismic rates. Nature, 427:436.
Di Toro, G., Han, R., Hirose, T., Paola, N. D., Nielsen, S., Mizoguchi, K., Ferri, F., Cocco, M.,
and Shimamoto, T. (2011). Fault lubrication during earthquakes. Nature, 471:494.
Doornhof, D., Kristiansen, T. G., Nagel, N. B., Pattilo, P. D., and Sayers, C. (2006). Compaction
and subsidence. Oileld Review, pages 50{68.
Dowell, N., Shah, P. S. F. N., and Maitland, G. C. (2017). The role of co
2
capture and utilization
in mitigating climate change. Nat. Clim. Change, 7:243249.
Elbanna, A. E. and Carlson, J. M. (2014). A two-scale model for sheared fault gouge: Compe-
tition between macroscopic disorder and local viscoplasticity. J. Geophys. Res. Solid Earth,
119:48414859.
Eliebid, M., Mahmoud, M., Shawabkeh, R., Elkatatny, S., and Hussein, I. (2017). Eect of co
2
content on the natural gas production from tight gas sandstone reservoirs. SPE Middle East
Oil & Gas Show and Conference, SPE.
Ergenzinger, C., Seifried, R., and Eberhard, P. (2011). A discrete element model to describe
failure of strong rock in uniaxial compression. Granul. Matter, 13:341364.
Eringen, A. (1994). A continuum theory of porous elastic solids. Int. J. Eng. Sci., 32:1337.
Espinoza, D., Vandamme, M., Dangla, P., Pereira, J. M., and Vidal-Gilbert, S. (2013). A trans-
verse isotropic model for microporous solids: Application to coal matrix adsorption and swelling.
J. Geophys. Res. Solid Earth, 118:6113{6123.
Ettles, C. M. (1986). The thermal control of friction at high sliding speeds. J. Tribol., 108:98104.
Ferdowsi, B., Gria, M., Guyer, R. A., Johnson, P. A., Marone, C., and Carmeliet, J. (2013).
Microslips as precursors of large slip events in the stick-slip dynamics of sheared granular layers:
A discrete element model analysis. Geophys. Res. Lett., 40:4194.
120
Finsy, V., Ma, L., Alaerts, L., DeVos, D. E., Baron, G. V., and Denayer, J. F. M. (2009).
Separation of co
2
/ch
4
mixtures with the mil-53(al) metal-organic framework. Microporous
Mesoporous Mater., 120:221{230.
Fossen, H., Schultz, R. A., Zhipton, Z. K., and Mair, K. (2007). Deformation bands in sandstone:
a review. J. Geol. Soc. London, 164:755{769.
Fraige, F. and Langston, P. (2004). Integration schemes and damping algorithms in distinct
element models. Adv. Powder Technol., 15:227246.
Frey, G. (2008). Hybrid porous solids: Past, present, future. Chem. Soc. Rev., 37:191{214.
Furukawa, H., Ko, N., Go, Y. B., Aratani, N., Choi, S. B., Choi, E., Yazaydin, A. O., Snurr,
R. Q.,
OKeee, M., and Yaghi, O. M. (2010). Ultra-high porosity in metal-organic frameworks.
Science, 329:424{428.
Gasparik, M., Bertier, P., Gensterblum, Y., Ghanizadeh, A., Krooss, B. M., and Littke, R. (2014).
Geological controls on the methane storage capacity in organic-rich shales. Int. J. Coal Geol.,
123:34{51.
Gelb, L., Gubbins, K. E., Radhakrishnan, R., and Sliwinska-Bartkowiak, M. (1999). Phase sepa-
ration in conned systems. Rep. Prog. Phys., 62:1573{1659.
Gensterblum, Y., Merkel, A., Busch, A., and Krooss, B. M. (2013). High-pressure ch
4
and co
2
sorption isotherms as a function of coal maturity and the in
uence of moisture. Int. J. Coal
Geol., 118:45.
Ghanbarian, B., Hunt, A., Ewing, R., and Sahimi, M. (2013a). Tortuosity in porous media: A
critical review. Soil Sci. Soc. Am. J., 77:1461.
Ghanbarian, B., Hunt, A., Sahimi, M., Ewing, R., and Skinner, T. (2013b). Percolation theory
generates a physically based description of tortuosity in saturated and unsaturated porous
media. Soil Sci. Soc. Am. J., 77:1920.
Ghassemzadeh, J., Hashemi, M., Sartor, L., and Sahimi, M. (2001). Pore network simulation
of
uid imbibition into paper during coating processes: I. model development. AIChE J.,
47:519535.
Ghassemzadeh, J. and Sahimi, M. (2004). Pore network simulation of
uid imbibition into paper
during coatingiii: Modeling of two-phase
ow. Chem. Eng. Sci., 59:2281.
Ghassemzadeh, J., Xu, L., Tsotsis, T. T., and Sahimi, M. (2000). Statistical mechanics and
molecular simulation of adsorption of gas mixtures in microporous materials: Pillared clays
and carbon molecular sieve membranes. J. Phys. Chem. B, 104:3892{3905.
Ghou, A. and Maurin, G. (2010). Hybrid monte carlo simulations combined with a phase mixture
model to predict the structural transitions of a porous metal-organic framework material upon
adsorption of guest molecules. J. Phys. Chem. C, 114:6496{6502.
Ghysels, A., Vanduyfhuys, L., Vandichel, M., Waroquier, M., Speybroeck, V. V., and Smit, B.
(2013). On the thermodynamics of framework breathing: A free energy model for gas adsorption
in mil-53. J. Phys. Chem. C, 117:11540{11554.
Giesting, P., Guggenheim, S., van Groos, A. F. K., and Busch, A. (2012). Interaction of car-
bon dioxide with na-exchanged montmorillonite at pressures to 640 bars: Implications for co
2
sequestration. Int. J. Greenhouse Gas Contr., 8:73{81.
121
Goldsby, D. and Tullis, T. (2002). Low frictional strength of quartz rocks at subseismic slip rates.
Geophys. Res. Lett., 29:1844.
G omez-Gualdr on, D., Col on, Y. J., Zhang, X., Wang, T. C., Chen, Y. S., Hupp, J. T., Yildirim, Y.,
Farha, O. K., Zhang, J., and Snurr, R. Q. (2016). Evaluating topologically diverse metal.organic
frameworks for cryo-adsorbed hydrogen storage. Energy Envi. Sci., 9:3279{3289.
Gor, G. and A. V. Neimark, A. (2011). Adsorption-induced deformation of mesoporous solids:
Macroscopic approach and density functional theory. Langmuir, 27:6926{6931.
Grosman, A. and Ortega, C. (2008). In
uence of elastic deformation of porous materials in
adsorption-desorption process: a thermodynamic approach. Phys. Rev. B, 78:085433.
Guo, L., Xiao, L., Shan, X., and Zhang, X. (2016). Modeling adsorption with lattice boltzmann
equation. Sci. Rep., 6:paper 27134.
Guo, Z., Shi, B., , and Wang, N. (2000). Lattice bgk model for incompressible navier-stokes
equation. J. Comput. Phys., 165:288.
Guyer, R. and Johnson, P. A. (2009). Nonlinear Mesoscopic Elasticity: The Complex Behavior
of Rocks, Soil, Concrete. Wiley.
Guyer, R. and Kim, H. A. (2015). Theoretical model for
uid-solid coupling in porous materials.
Phys. Rev. E, 91:042406.
Hassanizadeh, S. and Gray, W. (1979a). General conservation equations for multi-phase systems:
1. averaging procedure. Adv. Water Resour., 2:131.
Hassanizadeh, S. and Gray, W. (1979b). General conservation equations for multi-phase systems:
2. mass, momenta, energy, and entropy equations. Adv. Water Resour., 2:191.
Hassanizadeh, S. and Gray, W. (1990). Mechanics and thermodynamics of multiphase
ow in
porous media including interface boundaries. Adv. Water Resour., 13:169.
He, D., Ekere, N., , and Cai, L. (1999). Computer simulation of random packing of unequal
particles. Phys. Rev. E, 60:7098.
He, Y. B., Zhou, W., Qian, G. D., and Chen, B. L. (2014). Methane storage in metalorganic
frameworks. Chem. Soc. Rev., 43:56575678.
Heller, R. and Zoback, M. (2014). Adsorption of methane and carbon dioxide on gas shale and
pure mineral samples. J. Unconven. Oil Gas Resour., 8:14{24.
Hemmen, H. E. G. R., Fonseca, D. M., Hansen, E. L., Fossum, J. O., and Plivelic, T. S. (2012). X-
ray sudies of carbon dioxide intercalation in na-
uorohectorite clay at near-ambient conditions.
Langmuir, 28:1678{1682.
Herm, Z., Swisher, J. A., Smit, B., Krishna, R., and Long, J. R. (2011). Metal-organic frameworks
as adsorbents for hydrogen purication and precombustion carbon dioxide capture. J. Am.
Chem. Soc., 133:5664{5667.
Hirose, T. and Shimamoto, T. (2005). Growth of a molten zone as a mechanism of slip weakening
of simulated faults in gabbro during frictional melting. J. Geophys. Res., 110:B05202.
Hocker, T., Rajendran, A., and Mazzotti, M. (2003). Measuring and modeling supercritical
adsorption in porous solids. carbon dioxide on 13x zeolite and on silica gel. Langmuir, 19:1254{
1267.
122
Hol, S., Peach, C. J., and Spiers, C. J. (2011). Applied stress reduces the co
2
sorption capacity
of coal. Int. J. Coal Geol., 85:128{142.
Hol, S., Peach, C. J., and Spiers, C. J. (2012). Eect of 3-d stress state on adsorption of co
2
by
coal. Int. J. Coal Geol., 93:1{15.
Hol, S. and Spiers, C. J. (2012). Competition between adsorption-induced swelling and elastic
compression of coal at co
2
pressure up to 100 mpa. J. Mech. Phys. Solids, 60:1862{1882.
Hutzler, S., P eron, N., Weaire, D., and Drenckhan, W. (2004). The foam/emulsion analogy in
structure and rainage. Eur. Phys. J. E, 14:381.
Iliev, O., Mikeli c, A., and Popov, P. (2008). On upscaling certain
ows in deformable porous
media. Multiscale Model. Simul., 7:93.
Iritani, E., Katagiri, N., Yamaguchi, K., and Cho, J.-H. (2006). Compression-permeability prop-
erties of compressed bed of superabsorbent hydrogel particles. Drying Technol., 24:1243.
Jasinski, L., Sangar e, D., Adler, P., Mourzenko, V., Thovert, J.-F., Gland, N., and B ekri, S.
(2015). Transport properties of a bentheim sandstone under deformation. Phys. Rev. E,
91:013304.
Jeroy, M., Fuchs, A. H., and Boutin, A. (2008). Structural changes in nanoporous solids due to
uid adsorption: thermodynamic analysis and monte carlo simulations. Chem. Commun., page
32753277.
Jerauld, G., Hateld, J. C., Scriven, L. E., and Davis, H. T. (1984a). Percolation and conduction
on voronoi and triangular networks: a case study in topological disorder. J. Phys. C, 17:1519{
1526.
Jerauld, G., Scriven, L. E., and Davis, H. T. (1984b). Percolation and conduction on the 3d voronoi
and regular networks: a second case study in topological disorder. J. Phys. C, 17:3429{3436.
Ji, L., Zhang, T., and Zhang, K. L. M. J. Q. X. (2012). Experimental investigation of main
controls to methane adsorption in lay-rich rocks. Appl. Geochem., 27:2533{2545.
Jivkov, A., Hollis, C., Etiese, F., McDonald, S., and Withers, P. (2013). A novel architecture for
pore network modeling with applications to permeability of porous media. J. Hydrol., 486:246.
Johnson, K. (1985). Contact Mechanics. Cambridge University Press, Cambridge.
Kadoura, A., Nair, A. K. N., and Sun, S. (2016). Adsorption of carbon dioxide, methane and
their mixture by montmorillonite in the presence of water. Microporous Mesoporous Mater.,
225:331{341.
Kang, S., Fathi, E., Ambrose, R. J., Akkutlu, I. Y., and Sigal, R. F. (2011). Carbon dioxide
storage capacity of organic-rich shales. SPE J., 16:842{855.
Karada, E. (2010). Investigation of swelling/sorption characteristics of highly swollen aam/amps
hydrogels and semi ipns with peg as biopotential sorbent. J. Encapsulation Adsorption Sci.,
1:6.
Kerstein, A. (1983). Equivalence of the void percolation problem for overlapping spheres and a
network problem. J. Phys. A, 16:3071.
Kierlik, E., Monson, P. A., Rosinberg, M. L., Sarkisov, L., and Tarjus, G. (2001). Capillary
condensation in disordered porous materials: Hysteresis versus equilibrium behavior. Phys.
Rev. Lett., 87:055701.
123
Kierlik, E., Rosinberg, M. L., Tarjus, G., and Pitard, E. (1998). Mean-spherical approximation
for a lattice model of a
uid in a disordered matrix. Mol. Phys., 95:341{351.
Kim, H., Shi, Y., He, J., Lee, H. H., and Lee, C. H. (2011). Adsorption characteristics of co
2
and
ch
4
on dry and wet coal from subcritical to supercritical conditions. Chem. Eng. J., 171:45{53.
Kim, N., Harale, A., Tsotsis, T. T., and Sahimi, M. (2007). Atomistic simulation of nanoporous
layered double hydroxide materials and their properties: Ii adsorption and diusion. J. Chem.
Phys., 127:224701.
Kim, N., Kim, Y., Tsotsis, T. T., and Sahimi, M. (2005). Atomistic simulations of nanoporous
layered double hydroxide materials and their properties i. structural modeling. J. Chem. Phys.,
122:214713.
Kim, T., Cho, J., and Lee, K. S. (2017). Evaluation of co
2
injection in shale gas reservoirs with
multi-component transport and geomechanical eects. Appl. Energy, 190:1195.
Kitagawa, S. and Uemura, K. (2005). Dynamic porous properties of coordination polymers in-
spired by hydrogen bonds. Chem. Soc. Rev., 34:109{119.
Kitajima, H., Chester, F., and Chester, J. (2011). Dynamic weakening of gouge layers in high-
speed shear experiments: assessment of temperature-dependent friction, thermal pressurization,
and
ash heating. J. Geophys. Res.-Solid Earth, 116:B08309.
Kitajima, H., Chester, J., Chester, F., and Shimamoto, T. (2010). High-speed friction of dis-
aggregated ultra cataclasite in rotary shear: characterization of frictional heating, mechanical
behavior, and microstructure evolution. J. Geophys. Res., 115:408.
Knackstedt, M., Shappard, A. P., and Pinczewski, W. V. (1998). Simulation of mercury porosime-
try on correlated grids: evidence for extended correlated heterogeneity at the pore scale. Phys.
Rev. E, 58:R6923{R6926.
Knackstedt, M., Shappard, A. P., and Sahimi, M. (2001). Pore network modeling of two-phase
ow in porous rock: The eect of correlated heterogeneity. Adv. Water Resour., 24:257{278.
Koch, D. and Brady, J. (1985). Dispersion in xed beds. J. Fluid Mech., 154:399.
Koehler, S. A., Hilgenfeldt, S., and Stone, H. A. (2000). A generalized view of foam drainage:
Experiment and theory. Langmuir, 16:6327.
Kowalczyk, P., Furmaniak, S., Gauden, P. A., and Terzyk, A. P. (2010). Carbon dioxide
adsorption-induced deformation of microporous carbons. J. Phys. Chem. C, 114:5126{5133.
Kozeny, J. (1927). Ueber kapillare leitung des wassers im boden. Sitzungsber Akad. Wiss. Wien,
136:271.
Kruggel-Emden, H., Sturm, M., Wirtz, S., and Scherer, V. (2008). Selection of an appropriate time
integration scheme for the discrete element method (dem). Comput. Chem. Eng., 32:22632279.
Kulasinski, K., Guyer, R., Derome, D., and Carmeliet, J. (2015a). Poroelastic model for
adsorption-induced deformation of biopolymers obtained from molecular simulations. Phys.
Rev. E, 92:022605.
Kulasinski, K., Guyer, R., Derome, D., and Carmeliet, J. (2015b). Water adsorption in wood
microbril: Role of crystalline-amorphous interface. Biomacromolecules, 16:2972{2978.
124
Lachenbruch, A. (1980). Frictional heating,
uid pressure, and the resistance to fault motion. J.
Geophys. Res., 85:6097{6112.
Lai, W. and Mow, V. C. (1980). Drag-induced compression of tissues around a diarthrodial joint.
Biorheology, 17:111.
Lee, J., Weingarten, M., and Ge, S. M. (2016). Induced seismicity: the potential hazard from
shale gas development and co
2
geologic storage. Geosci. J., 20:137{148.
Leonowicz, M. and Oleszak, D. (2009). Kinetics of swelling and drug release from pni-
paam/aliginate stimuli responsive hydrogels. Solid State Phenom., 154:17.
Li, H., Eddaoudi, M.,
OKeee, M., and Yaghi, O. M. (1999). Design and synthesis of an excep-
tionally stable and highly porous metal-organic framework. Nature, 402:276{279.
Li, J. R., Sculley, J., and Zhou, H. C. (2012). Metalorganic frameworks for separations. Chem.
Rev., 112:869{932.
Lim, S., Ashby, M., and Brunton, J. (1989). The eects of sliding conditions on the dry friction
of metals. Acta Metallurgica, 37:767772.
Liu, H., Zhang, L., and Seaton, N. A. (1992). Determination of the connectivity of porous solids
from nitrogen sorption measurements. ii. generalisation. Chem. Eng. Sci., 47:4393{4404.
Liu, J., Peach, C. J., and Spiers, H. Z. C. J. (2015). Thermodynamic models for swelling of
unconned coal due to adsorption of mixed gases. Fuel, 157:151{161.
Liu, L. and Bhatia, S. K. (2013). Molecular simulation of co
2
adsorption in the presence of water
in single-walled carbon nanotubes. J. Phys. Chem. C, 117:13479{13491.
Liu, Y., Her, J. H., Dailly, A., Ramirez-Cuesta, A. J., Neumann, D. A., and Brown, C. M. (2008).
Reversible structural transition in mil-53 with large temperature hysteresis. J. Am. Chem. Soc.,
130:11813{11818.
Liu, Y. and Wilcox, J. (2012). Eects of surface heterogeneity on the adsorption of co
2
in micro-
porous carbons. Environ. Sci. Technol., 46:1940.
Lopez-Echeverry, J. S., Reif-Acherman, S., and Araujo-Lopez, E. (2017). Peng-robinson equation
of state: 40 years through cubics. Fluid Phase Equilibria, 447:39{71.
Lorenceau, E., Louvet, N., Rouyer, F., and Pitois, O. (2009). Permeability of aqueous foams.
Eur. Phys. J. E, 28:293.
Loring, J., Ilton, E. S., Chen, J., Thompson, C. J., Martin, P. F., B en ezeth, P., Rosso, K. M.,
and Schaef, A. R. F. H. T. (2014). In Situ study of co
2
and h
2
o partitioning between
na.montmorillonite and variably wet supercritical carbon dioxide. Langmuir, 30:6120{6128.
Loring, J., Schaef, H. T., Turcu, R. V. F., Thompson, C. J., Miller, Q. R., Martin, P. F., Hu, J.,
Hoyt, D. W., Qafoku, O., Ilton, E. S., Felmy, A. R., and Rosso, K. M. (2012). In situ molecular
spectroscopic evidence for co
2
intercalation into montmorillonite in supercritical carbon dioxide.
Langmuir, 28:7125{7128.
Lu, Y., Ao, X., Tang, J., Jia, Y., Zhang, X., and Chen, Y. (2016). Swelling of shale in supercritical
carbon dioxide. J. Natural Gas Sci. Eng., 30:268{275.
Lu, Z., Abdou, M., and Ying, A. (2001). 3d micromechanical modeling of packed beds. J. Nuclear
Mater., 299:101.
125
Lucas, A., Mangeney, A., and Ampuero, J. P. (2014). Frictional weakening in landslides on earth
and on other planetary bodies. Nat. Commun., 5:3417.
Lucier, A., Zoback, M., Gupta, N., and Ramakrishnan, T. S. (2006). Geomechanical aspects of
co
2
sequestration in a deep saline reservoir in the ohio river valley region. Environ. Geosci.,
13:85.
Lyu, Q., Ranjith, P. G., Long, X., Kang, Y., and Huang, M. (2015). A review of shale swelling
by water adsorption. J. Nat. Gas Sci. Eng., 27:1421.
Mair, K., Elphick, S. C., and Main, I. G. (2002). In
uence of conning pressure on the mechanical
and structural evolution of laboratory deformation band. Geophysical Research Letters, 29:1410.
Mase, C. and Smith, L. (1987). Eects of frictional heating on the thermal, hydrologic, and
mechanical response of a fault. J. Geophys. Res., 92:6249{6272.
Mason, J., Sumida, K., Herm, Z. R., Krishna, R., and Long, J. R. (2011). Evaluating metal-organic
frameworks for post-combustion carbon dioxide capture via temperature swing adsorption. En-
ergy Environ Sci., 4:3030{3040.
Mason, J. A., Veenstra, M., and Long, J. R. (2014). Evaluating metalorganic frameworks for
natural gas storage. Chem. Sci., 5:3251.
Masoud, H. and Alexeev, A. (2010). Permeability and diusion through mechanically deformed
random polymer networks. Macromolecules, 43:10117.
Masoudi, R. and Pillai, K. M. (2010). Darcys law-based model for wicking in paper-like swelling
porous media. AIChE J., 56:2257.
Maxwell, S., Campbell, J. S. E., and Quirk, D. (2008). Microseismic deformation rate monitoring.
SPE paper 116596.
Mazumder, S. and Wolf, K. (2008). Dierential swelling and permeability change of coal in
response to co
2
injection for ecbm. International Journal of Coal Geology, 74:123{138.
McDonald, T., Mason, J. A., Kong, X., Bloch, E. D., Gygi, D., Dani, A., Crocella, V., Giordanino,
F., Odoh, S. O., Drisdell, W. S., Vlaisavljevich, B., Dzubak, A. L., Poloni, R., Schnell, S. K.,
Planas, N., , Lee, K., Pascal, T., Wan, L. F., Prendergast, D., Neaton, J. B., Smit, B., Kortright,
J. B., Gagliardi, L., Bordiga, S., Reimer, J. A., and Long, J. R. (2015). Cooperative insertion
of co
2
in diamine-appended metal-organic frameworks. Nature, 519:303{307.
Meakin, P. (1998). Fractals, Scaling and Growth far from Equilibrium. Cambridge University
Press, London, 2nd edition.
Mellot-Draznieks, C., Serre, C., Surbl, S., and Ferey, G. (2005). Very large swelling in hybrid
frameworks: A combined computational and powder diraction study. J. Am. Chem. Soc.,
127:16273{16280.
Merey, S. and Sinayuc, C. (2016). Analysis of carbon dioxide sequestration in shale gas reservoirs
by using experimental adsorption data and adsorption models. J. Nat. Gas Sci. Eng., 36:1087.
Michels, L., Fossum, J. O., Rozynek, Z., Hemmen, H., Rustenberg, K., Sobas, P. A., Kalantzopou-
los, G. N., Knudsen, K. D., Janek, M., Plivelic, T. S., and daSilva, G. J. (2015). Intercalation
and retention of carbon dioxide in a smectite clay promoted by interlayer cations. Sci. Rep.,
5:paper 8775.
Mindlin, R. (1949). Compliance of elastic bodies in contact. J. Appl. Mech., 16:259.
126
Mindlin, R. and Deresiewicz, H. (1953). Elastic spheres in contact under varying oblique forces.
J. Appl. Mech., 20:327.
Mishra, P., Edubilli, S., Mandal, B., and Gumma, S. (2013). Adsorption of co
2
, co, ch
4
and n
2
on dabco based metal organic frameworks. Microporous Mesoporous Mater., 169:75{80.
Mishra, P., Uppara, H. P., Mandal, B., and Gumma, S. (2014). Adsorption and separation of
carbon dioxide using mil-53 (al) metal-organic framework. Ind. Eng. Chem. Res., 53:19747{
19753.
Molinari, A., Estrin, Y., and Mercier, S. (1999). Dependence of the coecient of friction on the
sliding conditions in the high velocity range. J. Tribol., 121:3541.
Monson, P. (2012). Understanding adsorption/desorption hysteresis for
uids in mesoporous
materials using simple molecular models and classical density functional theory. Microporous
Mesoporous Mater., 160:47{66.
Mota, J., Martins, D., Lopes, D., Catarino, I., and Bonfait, G. (2017). Structural transitions in
the mil-53(al) metal-organic framework upon cryogenic hydrogen adsorption. J. Phys. Chem.
C, 121:24252{24263.
Neimark, A., Coudert, F. X., Boutin, A., and Fuchs, A. H. (2010). Stress-based model for the
breathing of metal-organic frameworks. J. Phys. Chem. Lett., 1:445{449.
Neimark, A., Coudert, F. X., Triguero, C., Boutin, A., Fuchs, A. H., Beurroies, I., and Denoyel,
R. (2011). Structural transitions in mil-53 (cr): View from outside and inside. Langmuir,
27:4734{4741.
Neuzil, C. E. (2003). Hydromechanical coupling in geologic processes. Hydrogeol. J., 11:41{83.
O
Hara, K., Mizoguchi, K., Shimamoto, T., and Hower, J. (2006). Experimental frictional heating
of coal gouge at seismic slip rates: evidence for devolatilization and thermal pressurization of
gouge
uids. Tectonophysics, 424:109118.
Ottiger, S., Pini, R., Storti, G., and Mazzotti, M. (2008). Measuring and modeling the competitive
adsorption of co
2
, ch
4
, and n
2
on a dry coal. Langmuir, 24:9531{9540.
Ougier-Simonin, A. and Zhu, W. (2013). Eects of pore
uid pressure on slip behaviors: an
experimental study. Geophys Res Lett, 40:26192624.
Ozdemir, E. and Schroeder, K. (2009). Eect of moisture on adsorption isotherms and adsorption
capacities of co
2
on coals. Energy & Fuels, 23:2821{2831.
Pan, Z., Connell, L., and Camilleri, M. (2010). Laboratory characterization of coal reservoir
permeability for primary and enhanced coalbed methane recovery. International Journal of
Coal Geology, 82:252{261.
Pan, Z. and Connell, L. D. (2007). A theoretical model for gas adsorption-induced coal swelling.
Int. J. Coal Geol., 69:243{252.
Pan, Z. and Connell, L. D. (2011). Modelling of anisotropic coal swelling and its impact on
permeability behavior for primary and enhanced coalbed methane recovery. Int. J. Coal Geol.,
85:257{267.
Peng, D. Y. and Robinson, D. B. (1976). A new two-constant equation of state. Ind. Eng. Chem.
Fund., 15:59{64.
127
Perera, M. S. A., Ranjith, P. G., Choi, S. K., and Airey, D. (2011). The eects of sub-critical
and super-critical carbon dioxide adsorption-induced coal matrix swelling on the permeability
of naturally fractured black coal. Energy, 36:6442{6450.
Perry, J., Perman, J. A., and Zaworotko, M. J. (2009). Design and synthesis of metal-organic
frame-works using metal-organic polyhedra as supermolecular building blocks. Chem. Soc. Rev.,
38:1400{1417.
Pham, N., Voronov, R., Tummala, N., and Papavassiliou, D. (2014). Bulk stress distributions in
the pore space of sphere-packed beds under darcy
ow conditions. Phys. Rev. E, 89:033016.
Pini, R., Ottiger, S., Burlini, L., Storti, G., and Mazzotti, M. (2009). Role of adsorption and
swelling on the dynamics of gas injection in coal. J. Geophys. Res., 114:B04203.
Pitois, O., Lorenceau, E., Louvet, N., and Rouyer, F. (2009). Specic surface area model for foam
permeability. Langmuir, 25:97.
Porter, B., Zauel, R., Stockman, H., Guldberg, R., and Fyhrie, D. (2005). 3-d computational
modeling of media
ow through scaolds in a perfusion bioreactor. J. Biomech., 38:543.
Potyondy, D. and Cundall, P. (2004). A bonded-particle model for rock. Int. J. Rock Mech. Min.
Sci., 41:13291364.
Psarras, P., Holmes, R., Vishal, V., and Wilcox, J. (2017). Methane and co
2
adsorption capacities
of kerogen in the eagle ford shale from molecular simulation. Acc. Chem. Res., 50:11818{1828.
Qajar, A., Daigle, H., and Prodanovi c, M. (2016). The eects of pore geometry on adsorption
equilibrium in shale formations and coal-beds: lattice density functional theory study. Fuel,
163:205{213.
Qian, Y., d'Humieres, D., and Lallemand, P. (1992). Lattice bgk models for navier-stokes equation.
Europhys. Lett., 17:479.
Rao, Q. and Leng, Y. (2016). Molecular understanding of co
2
and h
2
o in montmorillonite clay
interlayer under co
2
geological sequestration conditions. J. Phys. Chem. C, 120:2642{2654.
Ravikovitch and Neimark, A. V. (2006). Density functional theory model of adsorption deforma-
tion. Langmuir, 22:10864{10868.
Ravikovitch, P., Bogan, B. W., and Neimark, A. V. (2005). Nitrogen and carbon dioxide adsorp-
tion by soils. Environ. Sci. Technol., 39:4990{4995.
Ravikovitch, P. I., Vishnyakov, A., and Neimark, A. V. (2001). Density functional theories
and molecular simulations of adsorption and phase transitions in nanopores. Phys. Rev. E,
64:011602.
Reichenbach, C., Kalies, G., Lincke, J., Lassig, D., Krautscheid, H., Moellmer, J., and Thommes,
M. (2011). Unusual adsorption behavior of a highly
exible copper-based mof. Microporous
Mesoporous Mater., 42:592{601.
Remy, B., Khinast, J., and Glasser, B. (2009). Discrete element simulation of free
owing grains
in a four-bladed mixer. AIChE Journal, 55:20352048.
Rice, J. R. (2006). Heating and weakening of faults during earthquake slip. J. Geophys. Res.,
111:B05311.
128
Rice, J. R. (2017). Heating, weakening and shear localization in earthquake rupture. Phil. Trans.
R. Soc. A, 375:20160015.
Rother, G., Ilton, E. S., Wallacher, D., Hau, T., Schaef, H. T., Qafoku, O., Rosso, K. M., Felmy,
A. R., Krukowski, E. G., Stack, A. G., Grimm, N., and Bodnar, R. J. (2013). Co
2
sorp-
tion to subsingle hydration layer montmorillonite clay studied by excess sorption and neutron
diraction measurements. Environ. Sci. Technol., 47:205{ 211.
Rothert, E. and Shapiro, S. A. (2003). Microseismic monitoring of borehole
uid injections: Data
modeling and inversion for hydraulic properties of rocks. Geophysics, 68:685{689.
Sahimi, M. (2003). Heterogeneous Materials II. Springer, New York.
Sahimi, M. (2011). Flow and Transport in Porous Media and Fractured Rock. Wiley-VCH,
Weinheim, 2nd edition.
Sahimi, M. and Arbabi, S. (1993). Mechanics of disordered solids. ii. percolation on elastic
networks with bond-bending forces. Phys. Rev. B, 47:703.
Sahimi, M. and Imdakm, A. O. (1991). Hydrodynamics of particulate motion in porous media.
Phys. Rev. Lett., 66:1169.
Saksala, T. (2016). Numerical study of the in
uence of hydrostatic and conning pressure on
percussive drilling of hard rock. Comput. Geotech., 76:120128.
Salimi, H., Pourjavadi, A., Seidi, F., EftekharJahromi, P., and Soleyman, R. (2010). New smart
carrageenan-based superabsorbent hydrogel hybrid: Investigation of swelling rate and environ-
mental responsiveness. J. Appl. Polymer Sci., 117:3228.
Savoji, M. T. and Pourjavadi, A. (2006). Partially hydrolized kappa carrageenan polyacrylonitrile
as a novel biopolymer based superabsorbent hydrogel: Synthesis, characterization, and swelling
behaviors. Polymer Eng. Sci., 46:1778.
Schre
er, B. (2002). Mechanics and thermodynamics of saturated/unsaturated porous materials
and quantitative solutions. Appl. Mech. Rev., 55:355.
Scuderi, M. M., Carpenter, B. M., Johnson, P. A., and Marone, C. (2015). Poromechanics of
stick-slip frictional sliding and strength recovery on tectonic faults. J. Geophys. Res. Solid
Earth, 120:68956912.
Seaton, N. A. (1991). Determination of the connectivity of porous solids from nitrogen sorption
measurements. Chem. Eng. Sci., 46:1895{1909.
Serre, C., Bourrelly, S., Vimont, A., Ramsahye, N. A., Maurin, G., Llewellyn, P. L., Daturi, M.,
Filinchuk, Y., Leynaud, O., Barnes, P., and F erey, G. (2007). An explanation for the very large
breathing eect of a metalorganic framework during co
2
adsorption. Adv. Mater., 19:2246.
Shaer, G. (2010). Long-term eectiveness and consequences of carbon dioxide sequestration.
Nat. Geosci., 3:464467.
Shukla, R., Ranjith, P., Haque, A., and Choi, X. A. (2010). Review of studies on co
2
sequestration
and caprock integrity. Fuel, 89:2651{2664.
Simon, C. M., Braun, E., Carraro, C., and Smit, B. (2017). Statistical mechanical model of
gas adsorption in porous crystals with dynamic moieties. Proc. Natl. Acad. Sci. U. S. A.,
114:E287{E296.
129
Sing, K., Everett, D. H., Haul, R. A. W., Moscou, L., Pierotti, R. A., Rouqu erol, J., and Siemieni-
wska, T. (1985). Reporting physisorption data for gas/solid systems with special reference to
the determination of surface area and porosity. Pure Appl. Chem., 57:603.
Siriwardane, H., Gondle, R., and Smith, D. (2009). Shrinkage and swelling of coal induced
by desorption and sorption of
uids: theoretical model and interpretation of a eld project.
International Journal of Coal Geology, 77:188{202.
Smeraglia, L., Billi, A., Carminati, E., Cavallo, A., Toro, G. D., Spagnuolo, E., and Zorzi, F.
(2017). Ultra-thin clay layers facilitate seismic slip in carbonate faults. Sci. Rep., 7:paper 664.
Sone, H. and Shimamoto, T. (2009). Frictional resistance of faults during accelerating and decel-
erating earthquake slip. Nat. Geosci., 2:705.
Spagnuolo, E., Plmper, O., Violay, M., Cavallo, A., and Toro, G. D. (2015). Fast-moving disloca-
tions trigger
ash weakening in carbonate-bearing faults during earthquakes. Sci. Rep., 5:paper
16112.
Stewart, M., Ward, A., and Rector, D. (2006). A study of pore geometry eects on anisotropy in
hydraulic permeability using the lattice-boltzmann method. Adv. Water Resour., 29:1328.
Succi, S. (2001). The lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford
University Press, London.
Sudibandriyo, M., Mohammad, S. A., Robinson, R. L., and Gasem, K. A. M. (2011). Ono-kondo
model for high-pressure mixed-gas adsorption on activated carbons and coals. Energy and Fuels,
25:3355{3367.
Suh, M. P., Park, H. J., Prasad, T. K., and Lim, D. (2012). Hydrogen storage in metal organic
frameworks. Chem. Rev., 112:782835.
Sumida, K., Rogow, D. L., Mason, J. A., McDonald, T. M., Bloch, E. D., Herm, Z. R., Bae,
T. H., and Long, J. R. (2012). Carbon dioxide capture in metal-organic frameworks. Chem.
Rev., 112:724{781.
Tafti, T. A., Sahimi, M., Aminzadeh, F., and Sammis, C. G. (2013). Use of microseismicity for
determining the structure of the fracture network of large-scale porous media. Phys. Rev. E,
87:032152.
Tahmasebi, P., Sahimi, M., Kohanpur, A., and Valocchi, A. (2016). Pore-scale simulation of
ow
of co
2
and brine in reconstructed and actual 3d rock cores. J. Pet. Sci. Eng., 155:21{33.
Talu, O. and Myers, A. L. (2001). Reference potentials for adsorption of helium, argon, methane,
and krypton in high-silica zeolites. Colloid Surf. A, 187188:83{93.
Tanaka, H., Hiraide, S., Kondo, A., and Miyahara, M. T. (2015). Modeling and visualization
of co
2
adsorption on elastic layer-structured metal-organic framework-11: Toward a better
understanding of gate adsorption behavior. J. Phys. Chem. C, 119:1535{1543.
Tang, X. and Ripepi, N. (2017). High pressure supercritical carbon dioxide adsorption in coal:
Adsorption model and thermodynamic characteristics. J. CO
2
Util., 18:189{197.
Tang, X., Ripepi, N., Stadie, N. P., Yu, L. J., and Hall, M. R. (2016). A dual-site langmuir
equation for accurate estimation of high pressure deep shale gas resources. Fuel, 185:10{17.
Thovert, J.-F. and Adler, P. (2011). Grain reconstruction of porous media: Application to a
bentheim sandstone. Phys. Rev. E, 83:056116.
130
Torquato, S. (2002). Random Heterogeneous Materials. Springer, New York.
Tsotsis, T., Patel, H., Naja, B. F., Racherla, D., Knackstedt, M. A., and m. Sahimi (2004).
Overview of laboratory and modeling studies of carbon dioxide sequestration in coal beds. Ind.
Eng. Chem. Res., 43:2887{2901.
Valiullin, R., Naumov, S., Galvosas, P., K arger, J., Woo, H. J., Porcheron, F., and Monson, P. A.
(2006). Exploration of molecular dynamics during transient sorption of
uids in mesoporous
materials. Nature, 443:965{968.
Vallejo, L. E. and Lobo-Guerrero, S. (2005). The elastic moduli of clays with dispersed oversized
particles. Eng. Geol., 78:163{171.
van Bergen, F., Spiers, C., Floor, G., and Bots, P. (2009). Strain development in unconned coals
exposed to co
2
, ch
4
and ar: eect of moisture. International Journal of Coal Geology, 77:43{53.
Vandamme, M., Brochard, L., Lecampion, B., and Coussy, O. J. (2010). Adsorption and strain:
the co
2
-induced swelling of coal. J. Mech. Phys. Solids, 58:14891505.
Vilarrasa, V., Makhnenko, R., and Gheibi, S. (2016). Geomechanical analysis of the in
uence of
co
2
injection location on fault stability. J. Rock Mech. Geotech. Eng., 8:805{818.
Watson, D. (1981). Computing then-dimensional delaunay tessellation with application to voronoi
polytopes. Comp. J., 24:167.
Weinstein, T., Bennethum, L., and Cushman, J. (2008). Two-scale, three-phase theory for swelling
drug delivery systems. part i: Constitutive theory. J. Pharmaceutical Sci., 97:1878.
Whitaker, S. (1999). The Method of Volume Averaging. Kluwer Academic, Dordrecht.
Wibberley, C. and Shimamoto, T. (2005). Earthquake slip weakening and asperities explained by
thermal pressurization. Nature, 436:689.
Wilmer, C., Leaf, M., Lee, C. Y., Farha, O. K., Hauser, B. G., Hippo, J. T., and Snurr, R. R.
(2011). Large-scale screening of hypothetical metal-organic frameworks. Nat. Chem., 4:83{89.
Woo, H.-J., Sarkisov, L., and Monson, P. A. (2001). Mean-eld theory of
uid adsorption in a
porous glass. Langmuir, 17:7472{7475.
Wu, H., Thibault, C. G., Wang, H., Cychosz, K. A., Thommes, M., and Li, J. (2016). Eect of
temperature on hydrogen and carbon dioxide adsorption hysteresis in an ultramiroporous mof.
Microporous Mesoporous Mater., 219:186{189.
Wu, Q., Andreopoulos, Y., Xanthos, S., and Weinbaum, S. (2005). Dynamic compression of highly
compressible porous media with application to snow compaction. J. Fluid Mech., 542:281.
Xu, L., Tsotsis, T. T., and Sahimi, M. (2001). Statistical mechanics and molecular simulation of
adsorption of ternary gas mixtures in nanoporous materials. J. Chem. Phys., 114:7196{7210.
Yaghi, O. M.,
OKeee, M., Ockwig, N. W., Chae, H. K., Eddaoudi, M., and Kim, J. (2003).
Reticular synthesis and the design of new materials. Nature, 423:705{714.
Yang, K., Lu, X., Lin, Y., and Neimark, A. V. (2011). Eects of co
2
adsorption on coal deformation
during geological sequestration. J. Geophys. Res., 116:B08212.
Yang, N., Liu, S., and Yang, X. (2015). Molecular simulation of preferential adsorption of co
2
over ch
4
in na montmorillonite clay material. Appl. Surf. Sci., 356:1262{1271.
131
Yang, Q., Liu, D., Zhong, C., and Li, J. R. (2013). Development of computational methodologies
for metal-organic frameworks and their application in gas separations. Chem. Rev., 113:8261{
8323.
Yang, X. and Zhang, C. (2005). Structure and diusion behavior of dense carbon dioxide
uid in
clay-like slit pores by molecular dynamics simulation. Chem. Phys. Lett., 407:427{432.
Yue, Y., Rabone, J. A., Liu, H., Mahurin, S. M., Li, M. R., Wang, H., Lu, Z., Chen, B., Wang, J.,
Fang, Y., and Dai, S. (2015). A
exible metal-organic framework: Guest molecules controlled
dynamic gas adsorption. J. Phys. Chem. C, 119:9442{4449.
Zang, J. and Wang, K. (2017). Gas sorption-induced coal swelling kinetics and its eects on coal
permeability evolution: model development and analysis. Fuel, 189:164{177.
Zhang, J., Clennell, M. B., Liu, K., Pervukhina, M., Chen, G., and Dewhurst, D. N. (2016).
Methane and carbon dioxide adsorption on illite. Energy and Fuels, 30:10643.
Zhou, W., Wu, H., Hartman, M. R., and Yildirim, T. (2007). Hydrogen and methane adsorption
in metal-organic frameworks: A high-pressure volumetric study. J. Phys. Chem. C, 111:16131{
16137.
Zhu, H., Dhall, A., Mukherjee, S., and Datta, A. (2010). A model for
ow and deformation in
unsaturated swelling porous media. Transp. Porous Media, 84:335.
Zhu, W. and Wong, T.-F. (1999). Network modeling of the evolution of permeability and dilatancy
in compact rock. Water Resour. Res., 104:2963.
Zoback, M. D. and Zinke, J. (2002). Production-induced normal faulting in the valhall and ekosk
oil elds. Pure Appl Geophys, 159:403{420.
132
Abstract (if available)
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation
PDF
Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data
PDF
Machine-learning approaches for modeling of complex materials and media
PDF
Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale
PDF
The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media
PDF
Methanol synthesis in a membrane reactor
PDF
Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating
PDF
Molecular- and continuum-scale simulation of single- and two-phase flow of CO₂ and H₂O in mixed-layer clays and a heterogeneous sandstone
PDF
Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection
PDF
Flow and thermal transport at porous interfaces
PDF
Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing
PDF
The study of CO₂ mass transfer in brine and in brine-saturated Mt. Simon sandstone and the CO₂/brine induced evolution of its transport and mechanical properties
PDF
Hydrogen storage in carbon and silicon carbide nanotubes
PDF
Oscillatory and streaming flow due to small-amplitude vibrations in spherical geometry
PDF
A flow-through membrane reactor for destruction of a chemical warfare simulant
PDF
Molecular dynamics studies of protein aggregation in unbounded and confined media
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
A hybrid adsorbent-membrane reactor (HAMR) system for hydrogen production
PDF
Impact of exposure to brine/CO₂ on the mechanical and transport properties of the Mt. Simon sandstone
PDF
A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs
Asset Metadata
Creator
Bakhshian, Sahar
(author)
Core Title
Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Chemical Engineering
Publication Date
02/12/2018
Defense Date
01/23/2018
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
adsorption,deformation,fluid flow,OAI-PMH Harvest,porous media
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Sahimi, Muhammad (
committee chair
), Nakano, Aiichiro (
committee member
), Tsotsis, Theodore T. (
committee member
)
Creator Email
bakhshia@usc.edu,bakhshian.s@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c40-473985
Unique identifier
UC11265760
Identifier
etd-BakhshianS-6026.pdf (filename),usctheses-c40-473985 (legacy record id)
Legacy Identifier
etd-BakhshianS-6026.pdf
Dmrecord
473985
Document Type
Dissertation
Rights
Bakhshian, Sahar
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
adsorption
deformation
fluid flow
porous media