3 Hydfrac Theses

INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typemiter face, while others may be from any type of cornputer printer. The quality of this reproduction is dependent upon the quality of the copy subrnitted. Broken or indistinct pflnt, cdored or poor quality illustrations

and photographs, print bfeedthrough, substandard margins, and improper alignment can advenely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Alsol if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, &arts)

are reproduced by

secüoning the original, beginning at the upper left-hand m e r and continuing

from left to right in equal sections with small overiaps. Photographs incîuded in the original manusui-pt have b e n reprodueed xerographically in this copy.

Higher quality 6" x 9" bbck and white

photographie pnnts are available for any photographs or illustrations appeanng

in this copy for an additional charge. Contact UMI directly to order.

Bell & Howell Information and Leaming 300 North Zeeb Road, Ann Arbor, MI 48106-1346 USA 8001521-0600

COMPUTATIONAL MODELLING OF

FRACTURE AND DAMAGE IN POROELASTIC MEDIA

by Abbas Tofangchi Mahy ari May, 1997

Civil Engineering and Applied Mechanics McGill University, Montreal, Quebec, H3A 2K6

A thesis submitted to the Faculty of Graduate Studies and Research

in partial fulfillment of the requirements for the degree of Doctor of Philosophy

Ocopy~ight 1997, A.T. Mahyari

I*I

National Library of Canada

Bibliothèque nationale du Canada

Acquisitions and Bibliographie Services

Acquisitions et sewices bibliographiques

395 Wellington Street Ottawa ON K1A ON4 Cartach

395, rue Wellington Ottawa ON K I A ON4 Canada

Vour fim Votm rei6mw

Our fi& N o m reWrmce

The author has granted a nonexclusive licence allowing the National Library of Canada to reproduce, loan, distribute or sell copies of this thesis in microform, paper or electronic formats.

L'auteur a accordé une licence non exciusive permettant à la Bibliothèque nationale du Canada de reproduire, prêter, distribuer ou vendre des copies de cette thèse sous la forme de tnicrofiche/film, de reproduction sur papier ou sur fonnat électronique.

The author retains ownership of the copyright in ths thesis. Neither the thesis nor substantial extracts fiom it may be p ~ t e or d otherwise reproduced without the author's permission.

L'auteur conserve la propriété du droit d'auteur qui protège cette thèse. Ni la thèse ni des extraits substantiels de celle-ci ne doivent être imprimés ou autrement reproduits saos son autorisation.

To my parents,

Abstract

The classical theory of poroelasticity focuses on the coupled response of fluid flow and elastic deformation of porous media saturated with either an incompressible or a compressible pore fluid. The Theory of poroelasticity has been successfully applied to examine time-dependent transient phenomena in a variety of natural and synthetic

materials, including geomaterials and biomaterials. The assurnption of elastic behaviour of the porous skeleton in these developments is a significant limitation in the application of this theory to brinle geomaterials which could exhibit non-elastic phenomena either in

the forrn of initiation and extension of discrete fractures, or in the form of initiation and evolution of continuum darnage in the porous skeleton. The computational methodology developed in this study examines the effect of development of such dofects (fracture or damage) on the fluid transport characteristics and the poroelastic behaviour of saturated

geomatenals. The finite element based computational models for fracture and damage phenornena examine two-dimensional plane strain and axisymmetric problems. The classical theory of linear elastic fracture mechanics is extended to examine the timedependent behaviour of local effects at the crack tip in poroelastic media. The numerical procedure accounts for the stress singularity of the effective stress field at the crack tip. The darnage model based on the concept of continuum damage mechanics, takes into

account the alteration of the stifiess and permeability characteristics of porous material due to developrnent of micromechanical damage in the porous skeleton. The isotropie darnage criteria governing the evolution of stifhess and permeabiiity parameters are

characterized by the dependency of damage parameters on the distortional strain invariant,

As the applications of the theory of poroelasticity diversifi, attention needs to be focused

on other aspects of importance. The class of transient and steady crack extension in poroelastic media is recognized as an area of interest in geomechanics applications and in

energy resources recovery fiom geological formations. A computational algorithm is developed to examine the transient quasi-static crack extension in poroelastic media where the temporal and spatial variations of boundary conditions governing the displacement, traction and pore pressure fields are taken into account in the incremental analysis. The path of crack extension is established by a mixed-mode crack extension criterion applicable to the porous fabric. The computational modelling of steady state crack extension in poroelastic media at constant velocity is also examined for the plane main problems. The finite element formulations of the goveming equations, which are velocity-dependent, are developed by employing the Galerkin technique. The poroelastic

behaviour of materiai depends on the propagation velocity at the crack tip. The computational schemes developed in this study followed an extensive procedure of verification via known analytical solutions to poroelasticity problems and for limiting cases of initial undrained (t +O+) and h a 1 drained (t ++m) elastic responses recovrred through analogous problems in classical elasticity.

Résumé La théorie classique de la poroélasticité porte surtout sur les effets entre le débit et la déformation élastique d'un milieu poreux saturé avec soit un liquide compressible ou noncompressible. La théorie de la poroélasticité a été appliquée avec succès pour examiner les effets transitoires dans plusieurs matériaux naturels et synthétiques, incluant les géornatériaw et les biomatériaux. Mais, lorsque l'on suppose que le skeleton poreux

possede un comportement élastique, celi impose une importante limitation dans l'application de cette théorie pour des géomatériaux fragiles puisque ces matériaux pourraient démontrer des caractéristiques nonélastiques soit par l'initiation ou la propagation de fissures ou par une initiation et une évolution de dommage dans le skeleton poreux.

La méthode développé a l'aide d'ordinateur pour cette recherche

examine les effets du développement de ces défauts (dommages ou fissures) sur les

caractéristiques du transport d'un liquide, ainsi que le comportement poroélastique dans les géomatériawc staturés. Des modèles basés sur les éléments finis pour l'analyse des fissures et des phénomènes de dommages servent à examiner l'élongation en deux dimensions ainsi que d'autres problèmes axisymétrique. La théorie classique de la mécanique de l'élasticité simple de fractures est prolongée pour examiner les effets locaux à la pointe de la fissure dans un milieu porew par rapport au temps. Pour évaluer le dommage, le modèle, basé sur le concept de la mécanique du dommage dans un continuum, prends en compte les changements dans la rigidité et la perméabilité du matériel poreux dû au développement de micro-dommage dans le skeleton poreux. Les critères de dommages isotropiques gouvernant l'évolution des paramètres de rigidités et de perméabilités sont caractérisés par une dépendance aux paramètres de dommages dans les déformations distortionnelles.

À mesure que les applications de la théorie de la poroélasticité se diversifient on doit

mettre le point sur d'autres aspects de plus grandes importances. La classe d'extension de fissures statiques et transitoires dans un milieu poroélastique est reconnue comme un domain important vis-a-vis l'application de la géornécanique et la récupération de ressources d'énergie provenant de formations géologiques. Un algorithme sur ordinateur a éte développé pour examiner l'extension quasi-statique de fissures dans un milieu poroélastique, où les variations spaciales et temporelle des conditions governant le déplacement, la traction et les champs de pression sont tenus en compte dans l'analyse.

Le chemin prise par la fissure est établis par un critère d'extension de la fissure en mode mixe applicable aux matériaux poreux. Un modèle à l'aide d'ordinateur de l'extension des fissures à une vitesse constante dans un milieu poroélastique est aussi examiné pour des problèmes à élongation simple. Les équations gouvernantes venant des éléments finis, étant dépendantes de la vitesse, sont formulées en employant la technique de

Galerkin.

Le Comportement poroélastique des matériaux dépend de la vitesse de

propagation de la pointe de la fissure. La méthode développée à l'aide d'ordinateur pour cette recherche a suivi un procédé de vérification en utilisant des solutions connues analytiques aux problèmes de poroélasticité ainsi que aux problèmes analogué dans la théorie d'élasticité simple pour les cas initiaux où le milieu élastique est non-égoutté (t +

)'0 et égoutté (t + +oo).

Acknowledgments

Foremost, the author praises God for giving him health, strength and patience during many years of study.

The author wishes to express his sincere gratitude to his thesis supervisor, Professor A.P.S. Selvadurai for his extensive guidance, encouragement and support during the course of this research. His invaluable comments and carehil review are greatly appreciated. Sincere thanks and appreciation also go to Dr. Keyvan Sepehr and Mr. Nick Vannelli.

M.Eng. for creating a fiiendly and stimulating environment in the Geomechanics Research Laboratory. Special thanks are due to Dr. R-V. Craster, Lecture- Department of Theoreticai Mechanics, University of Nottingham, England for his usefui discussions and advice resulting in the improvement of the work. The financial support provided by the Natural Sciences and Engineering Research Council of Canada through research gants awarded to Professor A.P.S. Selvadurai is aiso acknowledged.

Finally, the author is indebted to his family who always provided him with love, care and moral support during al1 stages of his study.

Table of Contents Abstract

iii

Résumé

v

Acknowledgments

vii

List of Figures

xii

List of Symbols

xv

List of Publications Resulting from the Thesis

K~X

INTRODUCTION AND LITERATURE REVIE;W

1

General

1

Theory of Poroelasticity

1

Applications of Poroefasticity

3

Non-Linear Brittie Behaviùur of Soi1 Skeleton

4

1.4.1 Linear Fracture Mechanks

7

1.4.2 Continuum Damage Mechanics

13

Objectives and Scope of the Research

16

Statement of Originality and Contributions

18

COMPUTATIONAL DEVELOPMENTS IN POROELASTIC MEDIA

20

Introduction

20

Classical Theories of Poroelasticity

21

Governing Equations of Poroelasticity

23

Computational Modelling of Poroelastic Media

27

2.4.1 Finite Element Method

27

2.4.2 Boundary Element Method

28

2.5

Finite Element Formulations 2.5.1 Galerkin Weighted Residual Method 2.5.2 Finite Element Approximations

2.5.3 Instantaneous Poroelastic Response 2.5.4 Finite Element Discretizations

2.5.5 Tirne Variation of Boundary Conditions 3.

VERIFICATION OF FINITE ELEMENT CODE

3.1

One-Dimensional Consolidation

3.2

Consolidation of a Poroelastic Sphere

3.3

Plane Strain and Axisymmetnc Consolidation

1.

FRACTURE MECHANICS OF POROELASTIC MEDIA

4.1

Introduction

4.2

Fracture Mechanics Concepts 4.2.1 Modelling of Crack Tip Behaviour

4.2.2 Evaluation of Stress Intensity Factors 4.3

Stationary Poroelastic Fracture 5.3.1 Verification Exercises

4.4

Indentation of a Cracked Poroelastic Half-Space 4.4.1 Cylindncal Crack 4.4.2 Penny-Shaped Crack

5.

TRANSIENT M O W G BOUNDARY PROBLEMS IN

POROELASTIC MEDIA 5.1

Introduction

5.2

Crack Extension Modelling

5.2.1 Critena for Onset of Crack Extension 5.2.2 Orientation of Crack Growth

5.2.3 Crack Modelling Techniques

5.3

A Computational Mode1 for Crack Extension

5.4

Intemai Indentation of a Penny-Shaped Crack

5.5

Indentation of a Poroelastic Half-Space

6.

STEADY MOVING BOUNDARY PROBLEMS IN POROELASTIC MEDIA

6.1

Introduction

96

6.2

Steady Crack Propagation

98

6.3

Goveming Equations

102

6.4

Finite Element Formulations

1 O3

6.5

Verification Exercises

1 O6

6.6

Wedging of a Poroelastic Crack

Ill

6.6.1 Plane Strain Moving Punch

6.6.2 Axisymmetric Penetration of a Rigid Shell 7.

DAMAGE MECHANICS OF POROELASTIC MEDIA

7.1

Introduction

7.2

P ~ c i p l e of s Damage Mechanics 7.2.1 The Darnage Variable

7.2.2The Net Stress Concept 7.3

Initiation and Evolution of Damage

7.4

Experimental Results with Indication of Damage

7.5

Muence of Damage on Poroelastic Parameters 7.5.1 Etastic Parameters 7.5.2 Penneability Characteristics

7.6

Finite Element Procedure

7.5

Indentation Problem of a Damageable Geomaterial

8.

CONCLUSIONS AND RECOMMENDATIONS

8.1

Summary and Conclusions

8.2

RecommendationsforFutureWork

REFERENCES APPENDIX A

Finite Eiement Formulation of Poroelastic Media

List of Figures Figure 1.1

Indentation of a brittle geomatenal

Figure 1.2

The stress field near the crack tip

Figure 2.1

Solid isoparametric element

Figure 2.2

Time variation of prescribed displacement and traction boundary conditions

Figure 3.1

Finite element discretization of one-dimensional consolidation

Figure 3.2

Pore pressure evolution at different depths

Figure 3.3

Time variation of surface settlernent

Figure 3.4

Finite element discretization of poroelastic sphere

Figure 3.5

Pore pressure evolution at the center of poroelastic sphere

Figure 3.6

Finite element discretization of the poroelastic half-space under plane strain or axisymmetric consolidation

Figure 3.7

Tirne-dependent pore pressure evolution

45

Figure 3.8

Degree of consolidation for proelastic layer

46

Figure 4.1

Fundamental modes of deformation at the crack tip region

51

Figure 4.2

Quarter point singular elements

55

Figure 4.3

Node arrangement for computation of the stress intensity factors

58

Figure 4.4

Finite element discretization of a semi-infinite crack

60

in an infinite poroelanic medium Figure 4.5

Time dependent behaviour of the flaw opening mode stress

62

intensity factor KI Figure 4.6

Surface indentation of a cracked poroelastic half-space

65

Figure 4.7

Finite element discretization of poroelastic ha1f-space

66

Figure 4.8

Effect of stress singularity and crack permeability on (a) consolidation 68 settlement of the indentor; and (b) pore pressure at the crack tip

Figure 4.9

Effect of geometrical characteristics of the cylindrical crack on the

69

consolidation settfement of the indentor Figure 4.10

Effect of geometrical characteristics of the penny-shaped crack on

the consolidation settlernent of the indentor

Figure 5.1

Idealization of cracks; (a) nodal separation using two or four coincident nodes; (b) smeared crack (ASCE, 1982)

Figure 5.2

Identification of crack orientation

Figure 5.3

Interaction conditions on the faces of the crack

Figure 5.4

Internai indentation of a penny-shaped crack

Figure 5.5

Finite eiement discretization of indentation problem

Figure 5.6

Variations of (a) the crack opening stress intensity factor; and (b) evolution of crack length for a prescnbed total load P

Figure 5.7

Time-dependent load relaxation and stress intensity factor for a stationary crack subjected to a prescribed displacement

Figure 5.8

Time-dependent load relaxation and time-dependent evolution of crack length for a prescribed displacement loading

Figure 5.9

Surface indentation of a poroelastic half-space

Figure 5.10

Finite element discretization

Figure 5.1 1

Crack extension patterns for two mesh discretizations A and B

Figure 5.12

Effect of crack extension on the degree of consolidation

Figure 5.13

Effect of undrained compressibility on the crack extension path

Figure 6.1

Steady crack growth problem examined by Yoffe (195 1)

Figure 6.2

Steady propagation of a crack in a poroelastic medium

Figure 6.3

Finite element discretization of steady crack extension

Figure 6.4

Crack propagation criterion for uniform loading of crack

Figure 6.5

Variation of the stress intensity factor KI

Figure 6.6

Pore pressure distribution ahead of crack tip for velocity al=l

Figure 6.7

Wedging of a poroelastic medium by ngid hder:;xs

Figure 6.8

Effect of crack geometry, in situ stresses, aiid propagltion velocity on the stress intensity factor for a point force wedging

Figure 6.9

Effect of crack and punch geometry on the stress intensity

xiii

71

factor KI for a wedging &id strip Figure 6.10

Rigid cylindrical shell penetrating through a poroelastic medium

Figure 6.1 1

Finite element discretization of the ngid shell penetratian

Figure 6.12

Variation of the stress intensity factor Kiwith radius a for a ngid smooth shell penetrating steadily through a poroelastic medium

Figure 7.1

Typical stress-strain behaviour and permeability evolution in brinle geomaterials

Figure 7.2

Representative element of virgin and damage state of material

Figure 7.3

Hypothesis of strain equivalence

Figure 7.4

The uniaxial compression behavior of concrete

(Spooner and Dougill, 1975) Figure 7.5

Variations of volumevic strain and permeability coefficient

126

as a fùnction of differential stress (Zoback and Byerlee, 1975)

Figure 7.6

Variations of darnage variable with material parameters

129

Figure 7.7

The stress-stmin behaviour; evolution of damage and permeability

133

characteristics in uniaxial compression of sandstone

(adopted €iomCheng and Dusseauit, 1993)

Figure 7.8

Finite element discretization of indentation problern

Figure 7.9

Degree of consolidation settlement of indentor

Figure 7.10

Extenr of damaged zone (with D > 0.05) in time

Figure 7.1 1

Pore pressure evolution at different depths

Figure 7.12

Evolution of damage variable at the edge of indentor

xiv

List of Symbols Latin Symbols

Symbol

Description

a

= Radius of circular disc or constant

QI

= a V K ; Normalized propagation velocity

A

=

Initial cross sectionai area

=Net cross sectional area = Radius of cylindrical crack = Boundary of dornain

= Skempton's pore pressure coefficient = Strain-Displacement matrix in

finite element formulation

= Transverse (shear) elastic wave velocity = Rayleigh wave velocity

= Generalized coeficient of consolidation = Coefficient of consolidation = Continuity up to the fint order differentiation = Coupling matrix related to interaction behveen soi1 and fiuid = Modified coupling matrix in steady problems

= Matrix associated with boundary conditions in steady problems = Damage variable = Critical damage variable

= Elasticity matrix of porous skeleton = Young's elastic modulus of porous fabric = Young's elastic modulus of damaged material = Compressibility matrix of fluid

= Modified compressibility matrix of fluid in steady problems

= Force vecton due to extemal tractions, body forces and flows = Fracture energy release rate = Critical fracture energy release rate = Permeability matrix

= Heaviside step hinction = First stress invariant = Coefficient of hydrauiic conductivity of porous medium = Coefficient of hydraulic conductivity of damaged matenal = stifiess matrix of the

soi1 skeleton

= Opening mode stress intensity factor

= Shearing mode stress intensity factor

= Cntical stress intensity factor (fracture toughness) = Crack length

= Slope of ramp fùnction for loading = unit out-normal

vector

= direction cosine of out-normal vector on boundary B = Interpolation function for the displacement field = Interpolation hinction

for the pore pressure field

= Total extemal load = Excess pore fluid pressure = Pore pressure at node K

= Vector of nodal pore pressure in finite elernent formulation

= Radial distance of location = Domain of region = Component of displacement vector

= Component of displacement vector in the iLhdirection at node J = time

= time increment

T

= Nonnalized tirne

V

= Propagation velocity in steady problems

w*

= Weighting hinctions in Galerkin technique

x, y, 2

= Cartesian coordinates

factor

Greek Symbols

Syrnbol

Description

a

= Biot's material parameter

P

= Biot's matenai parameter

Eij

= Component of strain tensor

Y

= Integration constant (>O and 4)

Yw

= Unit weight of pore fluid

K

= MY,

(6)

= Vector of nodal displacements of solid skeleton in finite element

formulation = Kronecker delta huiction (=l

if i=j; =O if i#j)

= Elastic shear modulus of porous skeleton = Elastic shear modulus of darnaged material

= Drained Poisson's ratio = Undrained Poisson's

ratio

= Specific discharge vector

= Angle of crack extension in local polar coordinates at crack tip = m a s density of matenal = Component of total stIess tensor = Component of total net stress tensor = Component of effective stress tensor

xvii

= Total in situ stress in geornatenals =Exted

normal traction or in situ stresses

= Local coordinate system for isoparametric elements

= Shear strain energy

= Volumetric strain in the compressible pore fluid = Displacement of indentor

List of Publications Resulting from the Thesis 1. Selvadurai A.P.S. and A.T. Mahyari, 1997, 'Computational Modeling of the

Indentation of a Cracked Poroelastic Half-Space', International Journcil of Fracture, (in press). 2. Selvadurai A.P.S. and A.T. Mahyari, 1997, 'Computational Modelling of Steady Crack Extension in Poroelastic Media ', International Journal of Soli& und Structures, (in review).

3. Mahyari A.T. and A.P.S. Selvadurai, 1997, 'Enhanced Consolidation in Brittle

Geomateriais Susceptible to Damage ,' International Journal of Mechanics of CohesiveFriction01 M d e rials (in press). 4. Mahyarî A.T. and A.P.S. Selvadurai, 1997, 'Interaction of a Sphencal Cavity and a

Circular Crack in a Poroelastic Medium', ldhCianadian Congress ofApplied k6echunics, CANCAM 97, June 1-6, Québec, 13- 14.

5. Mahyari A.T. and A.P.S. Selvadurai, 1997, 'Poroelastic Behaviour of a Rigid Anchor

Plate Ernbedded in a Cracked Geornaterial',

8International Symposium on Numerid

Models in Geomechanics, NUMOG 97, July 2-4, Montreal, 334-340.

6. Selvadurai A.P.S. and A.T. Mahyari, 1997, 'Some Moving Boundary Problems Associated with Poroelastic Media ', 41h International ConBrence on Moving Boundmies 97, Aupst 27-29, Ghent, Belgiurn, 216-224.

CHAPTER 1

INTRODUCTION AND LITERATClRE REVIEW

1.1 General

Fluid saturated geomaterials such as soils and rocks consist of a porous fabric of solid particles which form a deformable soil or rock skeleton the voids of which are filled with a pore fluid. The porous fabric can consist of an assemblage of individual particles which dcrives its structure and stiffness through intergranuiar contacts such as in sand or

alternatively, it can consist of a solid matenal which contains an interconnected network

of channels such as in rocks. The multiphase nature of saturated geomatenais makes their fundamental properties significantly different fiom other single phase materials (solid or

liquid). The mechanical behaviour of geomatenals is generally a function of the solid fabnc and the interaction arnong the solid and fluid phases.

1.2 Theory of Poroelrsticity

The first recognition of the importance of the multiphase nahue of geomaterials is generdly attributed to Terzaghi (1923). In the development of the theory of effective

stress, Terzaghi (1923) postuiated that when a saturated soil or rock is subjected to

extemal loads (e-g. loads induced by foundations), these loads are partly carried by the soil skeleton and partly by the pore fluid. The ability of the pore Buid to share the extemal loading is an important development in the undentanding of the mechanical behaviour of saturated geomaterials. The second important deveiopment is the innuence of the multiphase nature of the geomaterial on its time-dependent response to extemal loading. This second aspect is clearly demonstrated by Terzaghi (1923) in the development of classicai theory of soil consolidotion. Terzaghi (1923) postulated that when a low peneability soil is subjected to an extemai loading, this load can cause an increase in pore fluid pressure. The soi1 skeleton will initially undergo an immediate

deformation without the dissipation of the excess pore pressures. With tirne and depending upon the drainage conditions, the induced excess pore pressures uill dissipate. As the excess pore pressure dissipes. the loads are transferred to the soil skeleton to

maintain the system in equilibrium and to satisfy the mechanical and hydraulic compatibility. niis load transfer process wiil lead to M e r deformations of soi1 skeleton. Consequently, the consolidation process of saturated geomaterials is dominated by the time-dependent coupling between the deformation of the soil skeleton and the flow of pore fluid through the voids within the soil skeleton. The original developments of Tenaghi (1923) were restncted to the one-dimensional behaviour of saturated soils. Tenaghi (1923) assumed that the soi1 skeleton is isotropie and elastic, and both the pore fluid and the soil particles were incompressible. The flow

of water through the porous skeleton was described by Darcy's Law (1 856). Biot (194 1, 1955,

1956) extended this concept to include the three-dimensional effects,

compressibility of pore fluid and soil particles and anisotropic behaviour of the soil skeleton. The theories proposed by Terzaghi (1 923) and Biot (194 1, 1955, 1956) establish the basis for the chssical theory cfporoelasticity for a fluid saturated medium. The classical theory of poroelasticity foms a fiamework for the examination of engineering problems associated with fluid saniratecl elastic media. The govemiDg equations involve

the full coupling of theories of elasticity and fluid transport as applied to porous coatinua.

1.3 Applications of Poroelasticity

As documented by Selvadurai (1996), the classical theory of poroelasticity has been successfully applied to examine time-dependent transient phenomena in variety of natural and synthetic materials, including geomaterials and biomaterials. In geomechanics the theory of poroelasticity has been successfully applied to the examination of consolidation of soils and soft rocks. The analytical approaches to the application of Biot's theory of poroelasticity to the study of the consolidation behaviour of foundations, either on the surface or embedded within a fluid saturated medium, are given by McNamee and Gibson

(1960a, b), Agbezuge and Deresiewicz (1974), Chiarella and Booker (1975), Gaszynski and Szefer (1978), Selvadurai and Yue (i994), Yue and Selvadurai (1994, 1995) and Lan and Selvadurai (1996). In the context of geomechanics, problems associated with land

subsidence due to withdrawal of water or energy nsources such as oil and natural gas can also be approached via the theory of poroelasticity. In the treatrnent of such problems, computational approaches to poroelastic media have also been used quite successfully. Examples of these are given by Schrefler and Simoni (1987) and Lewis and Schrefler (1987).

In recent years, the theory of fluid saturated poroelastic media have gained attention in die context of thermally driven movement of fluids in saturated geological media. .b important consideration in these applications is the relative compressibility of the pore fluid in cornparison with the soi1 skeletal fabnc. Examples of such applications both conceniing analytical developments and computational modelling are given by Aboutit er al. (1985), Booker and Sawidou (1985) and Selvadurai and Nguyen (1995), Giraud

and Rousset (1996).

The classical theory of poroelasticity and its developments to include effects of large strain in the soi1 skeleton, irreversible deformation of the soi! skeleton and other tirnedependent phenomena have found applications in the study of many naturai and synthetic rnaterials. In particular, biological materials mch as bone, and other tissues such as

arteries, skin, etc. have ais0 k e n examined by considering such constitutive responses. A documentation of recent developments in these areas is given in the publication by Selvadurai (1 996). The use of purely mathematical methods for the development of analytical solutions for poroelasticity problems represents dificult exercises in particular due to the timedependency associated with the response of the geological medium. For this reason and in view of the interest in the application of the theory of poroelasticity to practical problems, attention has focused on the development of numencal methods, such as the finite element method and the boundary integral equation methods for the solution of problems

in poroelasticity. These numerical methods allow the development of approximate solutions even for complex geometries and boundq conditions and for nonlinear behaviour of the skeletal material. The finite element method has been the most wideiy used procedure in engineering app!ications. Sandhu and Wilson (1969) were the fint to apply finite element methods to the study of problems associated with consolidating geomaterials. Ghaboussi and Wilson (1973) and Booker and Small (1975) have developed f ~ t element e procedures for the analysis of problems associated with surface Ioading of semi-infinite media. Selvadurai and Gopal (1986j and Schrefler and Simoni (1 987) have used mapped infilnite elements to investigate the consolidation behaviour of saturated geornaterial regions of infinite extent. Lewis and Schrefler (1987) and Nguyen

and Selvadurai (1995) have used the finite element methods in connection with modelling of thermal consolidation in porous media. Applications of boundary element procedures to poroelasticity problerns are documented by Cleary (1977), Banerjee and Butterfield (1981), Kuroki and Onishi (1982), Cheng and Liggea (1984a), Brebbia et. a2 (1984), and

Dominguez (1 992).

1.4 Non-Linear Brittle Behaviour of Soi1 Skeleton

The assumption of elastic behaviour of the porous skeleton is a significant limitation in

the application of the classical theory of poroelasticity to geomaterials which could

exhibit irreversible and non-linear phenomena in the behaviour of the soil skeleton. With most naturally occurring bnttle geornatenals, these phenomena can range from the generation of microcracks (i.e. damage) to the development of macrocracks (i.e. fractures). The generation of these flaws in soil fabric can alter the deformability and permeability characteristics of the saturated geomatenals. As a result, the consolidation behaviour of bnnle saturated geornaterials such as rocks and overconsolidated clay can be influenced by the creation of such defects in the soi1 skeleton. The darnage process can result in the development of surface discontinuities in the fom of rnicrocracks andor volume discontinuities in the form of microvoids. At the scale of

microcracks, the damage phenomena results in a discontinuous medium. On the macroscale, however, damage can be modeled as variables applicable to a continuum region (Kachanov, 1958). In contrast to continuum darnage phenomena, the fracture

process is localized at the crack tip and gives rise to discontinuous fields for the displacement, traction and pore pressure variables.

In generd, most of bnttle geomaterials experience simultaneously both foms of continuum damage and discrete fracture processes. The existence of microcracks in the

region near the crack tip so called process zone (Bazant, 1991) with a non-linear material

behaviour is an indication of this phenomena. Factors ùitluencing the dominant mode of flaw generation should include the effects of the state of stress, rate of loading,

microstructure of the geomaterial, presence of stress singularities (e.g. sharp contacts) and ability of a flaw to open and close. Universal criteria cannot be postulated to determine which mode of flaw generation governs the nodinear elastic behaviour of brittle

geomaterials. Experimental procedures should also be employed to examine the response of materials. However, the notion of continuum darnage cm be more relevant to semi-

brinle geomaterials such as soft rocks, overconsolidated clay and other porous geological media where pre-peak progressive softening can occw due to generation of rnicrovoids or microcracks. The phenomena of discrete cracking is expected to be more relevant to

discrete r

Figure 1.1 Indentation of a brinle geomaterial. predominantly brittle behaviour of geornatexials such as cornpetent rocks (e.g. granite) under conditions of low confining stresses.

Classicd linear elastic fracture mechanics (see e.g. Broek, 1982; Bazant, 1991) cari be adopted in fracture mechanics analysis associated with poroelastic media. The initiation and extension of discrete hctures is govemed by the local effective stress fields near the

crack tip in poroelastic media The extension of a crack in a fluid sanirated material gives

rise to both transient and steady movhg boundary problems in poroelastic media. in these problems, the boundary conditions governing the displacements, tractions and pore fluid pressures change spatially and temporally as the crack extends (Figure 1.1). The pore pressure boundary conditions or drainage conditions on the faces of the crack cm be specified for different physical situations corresponding to permeable and impermeable

fracture surfaces.

The development of microcracks and microvoids in the porous skeleton of saturated geomaterials has an immediate eEect altering the elastic stiffiiess of the porous skeleton.

The secondq influence of nich damage processes manifests in the fonn of alteration of

the permeability characteristics of the porous medium (Figure 1.1). As a result, the fluid migration characteristics and the rate of consolidation within the saturated geomaterial

can be altered. In contrast to discrete fractures, the micromechanical damage of porous skeleton does not result in any discontinuity in displacement, traction and pore pressure fields within the porous medium. Also the darnage effects are govemed by the global state of stresses in the porous medium. The effect of soi1 skeletai damage on the timedependent poroelastic behaviour of saturated geomaterials can be examined by incorporating the concept of continuum darnage mechanics (Kachanov, 1958) into the classical theory of poroelasticity. It can be achieved by representing the stiffness properties and permeability characteristics of a porous medium as a fûnction of the state of darnage in the material.

1.4.1 Linear Fracture Mechanics

The linear elastic fracture mechanics (LEFM) is based on the theories of fracture originaily proposed by Griffith (192 1, 1924). He examined the stress field near elliptical flaws and postulated a critenon which defines the conditions necessary for the extension

of a crack. He developed a fracture theory based on the concept of energy balance and proposed that a flaw becornes unstable when the elastic strain energy change which results fiom an increment of crack growth, is larger than the surface energy (i.e. resistance) of material. The assurnption of linear elastic behaviour is central to the developrnents proposed by Griffith (1921, 1924). The Griffith concepts were extended by Irwin (1957, 1958) to examine the behaviour of the local effects at the tip of shaq cracks.

lrwin (195 7, 1958) adopted the analytical technique developed by Westergaard (1 938) in analysis of stresses and displacements ahead of a sharp crack in elastic materials and showed that the crack tip solutions couid be described by a single constant. This constant

which characterises the displacement and stress fields near the crack tip, later became

known as the stress intensity factor.

crack

-

crack tip

Figure 1.2 The stress field near the crack tip.

The assumption of linear elasticity in elastic fracture mechanics allows the identification of various modes of deformation of the crack tip region which can initiate the growth of cracks in poroelastic materials under the influence of effective stress fields. Invin (1957, 1958) observed that in general there are three independent modes of deformation at the

crack tip region for eiastic materials. Three modes of fiachue can in general be classified as mode 1 (in-plane opening of crack), mode II (in-plane shearing of crack) and mode III

(anti-plane shearing of crack). The influence of extemal stress States cm be represented as a linear combination of these three modes of deformation at the crack tip.

Inÿin (1957, 1958) showed that under plane strain conditions near the crack tip, asymptotic stress field a, (Figure 1.2) in an isotropie elastic material takes the following

form:

where

0,

is the stress tensor, r and 0 are the polar coordinates at the reference system

l;, is a dirnensionless function of 0 which depends only on the geomeûy of the crack. The stress intensity factor K charactenses the crack tip located at the crack tip (Figure 1.2).

behaviour. Equation 1.1 indicates a sfress singuIarity of order r"R at the crack tip as r+O. It is shown that the displacement field near the crack tip varies with an order of r".

In poroelastic materids, the behaviour of fractures is governed by the time-dependent pore pressure diffusion process at the crack tip region which depends on the pore pressure boundary conditions on the crack faces. As the pore pressure drainage takes place. the sness intensity factors at the crack tip Vary fiom the limiting case of elastic undrained behaviour (with Poisson's ratio vu) to the elasticfuly drained behaviour (with Poisson's ratio v ) . This tirne-dependent behaviour of the effective stress state in cracks is a characteristic feature of fracture mechanics problems associated with poroelastic media. Simons (1977) has shown that the stress singularity of order r-ln is preserved for the effective stress field at the crack tip in poroelastic media. This stress singularity is maintained throughout the time-dependent poroelastic behaviour of sahirated materiais. Craster and Atkinson (1 99 1) have shown that the pore pressure behaviour at the crack tip

as r+ O is spatially nonsingular for poroelastic fracture problems. The pore pressure field exhibits temporal singularity only for undrained behaviour of poroelastic material corresponding to time PO*, and as the pore Buid pressure dissipates at the crack tip. the pore pressure field becomes regular. They have s h o w that the pore pressure grudiena

at

the crack tip are, however, singular for the permeable pore pressure boundary conditions

on the crack faces. Applications of the theory of poroelasticity have largely focused on problems related to conventional geomechanicd applications where the emphasis has been on the evaluation

of deformations and consolidation rates in such media subjected primarily to extemai loads. As the applications of the theory diversify, attention needs to be focused on other

aspects of importance. The initiation and extension of cracks in poroelastic media are recognized as an area of both practical and fundamental interest particularly to endeavors involving hydraulic hchiring and energy resources recovery. The crack problems associated with fiuid saturatecf materials can be categorized into two basic classes: transient and steady state crack problems. In transient crack problems where the crack can be either stationary with fixed boundary conditions or can extend with moving boundary

conditions, the behaviour of a crack is tirne-dependent. With steadily propagating cracks. the poroelastic behaviour of the crack c m be examined by considering a "tirne independent" analysis where the time dependency is removed by appeal to a transformation which contains the velocity of steady crack extension.

i. Transient Cracks A great majority of applications of poroelasticity focus purely on initial boundary value

problerns where the boundary conditions are kept fixed both spatially and temporally.

The boundary conditions can relate to either the tractions or displacements exerted on the porous skeleton, or conditions prescribed on the pore fluid pressure and its spatial derivatives. The extension of cracks in poroelastic media generally results in a mixedmode moving boundary problem where the boundary conditions themselves are timedependent.

The class of transient fracture problems where there is the interaction of pore fluid flow and soi1 skeletal deformation, are important to enhanced energy resources recovery techniques which involve controlled hydraulic fiacturhg (Boone and Ingraffea, 1988; Detournay and Cheng, 1991). In hydraulic fiacturing techniques, during the initiation of fracture by the pressuization of the fluid, the pore Buid tends to leak out into the porous

formation and induce a volume expansion of material around the crack tip leading to partial closure of the fracture. Field evidence (Smith, 1985; and Nierode, 1985) suggests that the injection pressure required to initiate the fracture in poroelastic media is higher

than those predicted by classical elasticity models which neglect the influence of poroelasticity.

ii. Steady State Cruch

The problem of steady state constant velocity crack growth in an isotropic elastic medium was first examined by Yoffe (195 1) for the case of plane main deformation. Yoffe (195 1)

developed complete anaiflical resuits (i.e. solutions for displacement and stress fields) for a crack of fixed length propagating in an elastic medium, ahich is subjected to uniform remote tensile stresses, by a solution of the governing elastodynamic equations. Radok

(1956) and Broberg (1960) extended this study to examine the steady self-similar extension of a crack in an elastic material. The lirnits of propagation velocity in these studies were established in relation to Rayleigh wave velocity in elastic materiais. A recent review of these developments is presented by Freund (1990). The classicd fracture mechanics considerations for elastic materials indicate that when

the energy release rate of the fracture is greater than that of the driving mechanism, the extension of the crack is unstable and the failure is rapid. In poroelastic materials and other geomaterials which exhibit dissipative phenomena, the propagation of a fracture is likely to be dynamic and unsteady in the initial stages. However, a state of steady crack extension can be obtained at limiting times when the steady crack extension process has penisted over a long penod of t h e . For poroelastic materials, the flow of energy into the pore fiuid tends to stabilize the crack growth and that can result in a quasi-static extension of the crack with a certain velocity. The development of landslides in overconsolidated clay (Palmer and Rice, 1973; and Rice and Cleary, 1976) and the aftenhock events in an earthquake (Booker, 1974) have been attributed to processes associated with quasi-static crack propagation. The quasi-static crack growth also govems the hydraulic hcturing phenomena used quite extensively in oil resources recovery (Ruina, 1978; Huang and Russell, 1985a, b; Boone and Ingraffea, 1988). While the class of two-dimensional plane

strain problems give rise to steady state crack extension phenornena, the equivalent class of problems involving extension of circular cracks do not yield a steady state. Both the mathematical and computational modelling procedures have been applied to the solution of fracture mechanics problems in poroelastic media. Analytical approaches

provide the solution for simple fracture problems where the boundary conditions are kept fixed in time and space. Numerical procedures however allow the examination of

problems where the geometries are complex and the boundary conditions can Vary in space and time.

The problem of steady crack propagation in poroelastic media has been examined by Rice and Cleary (1976) and Rice and Simons (1976). They have found that consideration of

coupling between pore fluid flow and soi1 skeletal deformation c m result in some

poroelastic effects which prevent opening of cracks. The problem of the steady fracture in a poroelastic material was investigated by Rudnicki (1985) in c o ~ e c t i o nwith fault

initiation and propagation in porous rocks. These problems were examined in view of

their potential use in the study of earthquake mechanisms at the source location. Vandamme et al. (1989) used the displacemenr discontinuity method to examine the

transient plane strain behaviour of a stationary poroelastic fracture subjected to a sudden constant pore fluid pressure. They decornposed the fluid pressure loading into prescribed pore pressure and normal traction boundary conditions on the crack faces and examined

the time-dependent behaviour of the crack separateiy under each class of boundary conditions. Detoumay and Cheng (1991) have examined a similar problem in greater details to evaluate the design parameters, such as the time-dependent variation of the crack opening stress intensity factor, which are important to controlled processes of

hydradic hctwing.

The analytical solutions to problems of a semi-infinite crack propagating quasi-statica1l:-

in shear mode through a porous medium for various pore pressure boundary conditions on the crack faces have k e n given by Rice and Simons (1976) and Simons (1977). Cheng and Liggett (1984b) applied the boundary integral equation method to solve a

similar problem. Ruina (1978) examined the problem of a hydraulically loaded fkacture propagating at a constant rate with impermeable crack faces. He has s h o w that the apparent matenal resistance to fracture propagation, or apparent fracture toughness. increases with fracture propagation velocity . Huang and Russell ( 198ja, b) extended

these results to hydmulically loaded fractures propagating at a constant velocity where the fracture surfaces are permeable. Recently Atkùison and Craster (1991) and Craster and Atkinson (1996) have examined the problems of stationary and steadily propagating semi-infinite cracks embedded in onginally intact poroelastic media. These studies include the examination of crack behaviour in the presence of variable pore pressure boundary conditions at the crack faces and their potential influence on the pore pressure field and stress intensity factors at the

crack tip.

The numerical treatment of cracks occwing in poroelastic geomaterials bas received only limited attention. The computational modelling of fiacture problems in poroelastic media can be approached by adopting either finite element techniques or boundary element

techniques. Henshell and Shaw (1975) and Barsoum (1976) independently introduced the quarter point singularity elements in fracture analysis of elastic materials. Barsoum (1976) has s h o m that the i l R stress singularity can be incorporated into the quadratic isoparametric elements by shifiing the mid-side node close to the crack tip to their quarter point. The quarter point crack tip element offers its preference over any other special crack tip elements (such as hybrid elements etc., Owen and Fawkes, 1983) due to its simplicity and its ready availability in any finite element code. This element has been

successfully utilized to model crack problems in classical elasticity via finite element methods (IngrafYea, 1977a, b; Owen and Fawkes, 1983; Murti and Valliappan. 1986) and boundary element methods (Cruse and Wilson, 1977; Blandford et al., 198 1 ; Smith and Mason, 1982; Selvadurai and Au, 1989; Selvadurai and ten Busschen, 1995).

In order to develop the computational modeiling of the quasi-static crack extension in porous media, it is necessary to establish a crack extension critenon applicable to the solid skeleton. The subject of fracture extension in brinle elastic solids has been studied very extensively in recent years. Such studies have been motivated by the interest in the examination of crack extension in geomaterials such as concrete, rock and ice. Extensive accounts of these developments can be found in the literature on fracture mechanics (see. e.g. Liebowitz, 1968; Atkinson, 1979; Broek, 1982; Sih, 1991). In studies related to the

crack extension in brittle elastic materials, it is necessary to postulate two criteria The crack extension criterion establishes the stress conditions necessary for the onset of crack extension. The second relates to the criterion which establishes the orientation of crack growth. Boone and lngraffea (1988, 1990) have exarnined the two-dimensional extension of a plane crack in a poroelastic material which is dnven by the pore fiuid pressures. They developed a crack extension model based on the separation of the interface elements placed dong o priori known crack path when the maximum effective tensile stresses at the crack front reach a criticai value. This procedure neglects the fracture resistance of material which can be significant in fiacture phenomena of poroelastic media in the absence of in situ stresses.

1.4.2 Continuum Damage Mechanies

The concept of continuum damage mechanics (CDM)was developed by Kachanov (1958) to examine the tertiary creep of solids. This theory examines the effect of

microcrack development or any other microdefects on the behaviou. of materials pnor to

the development of macrocracks (Le. fractures). The process of damage is expected to be highly anisotropic in nature and could be restricted to localized zones. The theory of continuum damage mechanics has been widely used to predict the nonlinear response of a variety of materials including metals, concrete, composites, ice,

and geological matenals (Krajcinovic and Fonseka, 1981; Simo and Ju. 1987; Selvadurai and Au, 1991;Cheng and Dusseault, 1993; Hu and Selvadurai, 1995: etc.). The nonlinear behaviour of materials is modelled by introducing local damage variabies in the anaiysis.

Damage variables reflect average material degradation at macro-scale normally

associated with the classical continuum description. This facilitates the adaptation of the damage concept in the theov of elastiuity or in any other theory associated with classical continuum mechanics. The coupling of elasticity and damage models has been investigated by Sidoroff (MO), Krajcinovic (1984), Chow and Wang (1987). and

Lemaitre and Chaboche (1 990). The occurrence of microcracks and microvoids in the porous skeleton will have the immediate eEect altering the elastic stiffhess of the porous sanirated geomaterials. The secondary influence of such damage will manifest in the fom of alteration of the permeability charactenstics of the porous medium. The graduai degradation in the constitutive properties and evolution of hydraulic conductivity of the material are as a result of continuing growth of either dready existing microdefects or the progressive nucleation of new microdefects. For a given state of stress, the extent of darnage is an intdnsic property of the material which is defmed by the "damage evolution law". The criteria goveming the evolution of elasticity and permeability parameters as a result of damage can either be postuiated by appeal to micromechanics or detennined by

experiments.

Sima and l u (1987) have successhilly applied the concept of damage mechanics to simulate the experimentai r e d t s on concrete obtained by Scavuao et al. (1983). They developed an isotropie elasto-plastic damage mode1 employing a tensorid form of

damage variables. Cheng and Dusseault (1993) used an anisotropic darnage rnodel, using a vectorial representation of damage variable based on experimentai observations, to model the stress-strain behaviour of rocks and concrete. They postulated a darnage evolution law as a function of shear strain energy and exarnined the consolidation behaviour of a strip foundation on a poroelastic half-space which experiences the soil skeletal damage with no alteration in penneability characteristics. The development of damage, including initiation and codescence of microcracks, gives rise to nonlinear phenomena in the constitutive behaviour of fluid saturated geomaterials.

The effect of damage on either the degradation of elastic moduli or strength (in the form

of strain softening) of geomaterials such as rocks has been observed by Cook (1965) and Bieniawski e t al. (1967). Spooner and Dougill(1975) have observed similar phenomena on experîments conducted on concrete. The effect of microcrack generation on the

pemeability characteristics of sanirated geomaterials such as rocks and concrete has also been observed by Zoback and Byerlee (1 975), Samaha and Hover (1992), Shiping et al. ( 1994), and Kiyama er

al. (1 996).

1.5 Objectives and Scope of the Research

Saturated geological rnaterials with brinle behaviour can experience non-elastic phenomena in the soil skeleton either in the form of continuum micromechanical damage by development of rnicrocracks, or in the form of discrete fractures. Development of such

defects in the porous skeleton of saturated geomaterials will influence the fluid transport characteristics and the poroelastic behaviour of such materials. One of the objectives of

this research is to develop a methodology for computatioaal modelling of fracture and damage phenomena in porous media saturated with compressible pore fluids for both two-dimensional plane strain and axisymmetric situations. A great rnajority of applications of the theory of poroelasticity f o c w s purely on initial

boundary value problems where the boundary conditions are kept fixed both spatially and

temporally. The boundary conditions can relate to either the tractions exerted by the porous skeleton, or the displacements of the porous elastic skeleton or the pore fhid transport. There are a significant number of problems which specifically relate to rnoving boundmïes associated with hcture phenornena where the boundary conditions themselves are tirne dependent. 'Ihe class of time-dependent moving boundary problems associated nith either the transient quasi-static crack extension or the steady state crack extension in poroelastic media is exarnined in this thesis.

It is clear from the literature review that computational treatment of both non-elastic brinle behaviour of the porous fabric (either in forrn of continuum damage or in form of development of discrete cracks) and moving boundary problems arising from the transient and steady state extension of cracks in saturated geomaterials have received limited attention. Specifically, attention will be focused on the following aspects of the modelling and computational developments: i) Extension of the linear elastic fracture mechanics (LEFM)into poroelastic media. This can be achieved by the incorporation of a r'IR singularity in the effective stress field at the

crack tip in numencd modelling and the development of a computational procedure for the evaluation of time-dependent variation of the stress intensity factors at the crack tip located in poroelastic media. The verification of the computational procedures can be achieved by cornparison of the numerical results with analytical solutions availabie in the literature for simple geometries and boundary conditions.

ii) Development of a computational algorithm for tirne-dependent quasi-static crack extension in poroelastic media. The displacement, traction and pore water pressure

boundary conditions are updated as crack extension takes place through the discretized region. The path of crack extension c m generally be established by mixed-mode crack extension criteria applicable to porous fabrics.

iii) Examination of moving boundary problems in poroelasticity where the temporal and spatial variations of prescribed or fiee boundary conditions governing displacement, traction and pore pressure variables can be coasidered. The accuracy of the computational scheme can be verified with known analytical solutions for limiting elastic responses of initial undrained (r +O+)

and final drained (r ++a) solutions recovered through

analogous problems in classical elasticity. iv) Development of a computational formulation for the steady state constant velocity

crack extension in poroelastic media under conditions of plane strain and axial symmetry. The conventionai equations governing poroelasticity are modified to replace the tirne variable with a modified time variable which incorporates the constant velocity of the moving crack tip. The computationai procedures cm be verified by appeai to analytical solutions to the displacement and pore pressure fields at the crack tip. The steady moving boundary problems associated with wedging of cracks in saturated porous media can be

exarnined.

V)

Development of an isotropic darnage model which accounts for the alteration in

constitutive properties of the soi1 skeleton and permeability characteristics of the porous medium. The isotropic damage evolution criteria goveming the evolution of elastic stifhess and permeability parameters can be characterized by the dependency of damage variables on the distortional strain invariant. The damage model cm be applied to

examine the extent to which the poroelastic behaviour of saturated geomatends cm be infiuenced by the evolution of darnage in the porous skeleton. 1.6 Statement of Originality and Contributions 1. The work presented in this thesis extends the developments into the computational modelling of poroelastic media, to include instances where cracks and other defects in the

form of darnage evolves during the application of loading. To the authors knowledge. these extensions are considered to be novel and specific problems examined are original.

2. The work is extended to include cornputational rnodelling of steady state and transient effects associated with extension of cracks in fiuid saturated geomaterials. The work also includes new developments in moving boundary problems involving steady penetration of inclusions in poroelastic media.

3. n i e work described in the thesis has significant potential applications in geomechanics, geotechnical engineering, materials engineering, and energy resources recovery fiom geological formations. 1. The contributions resulting fiom the thesis have been published or accepted for

publication in leading International Journals and refereed conference proceedings with a

high degree of selectivity and standards.

CHAPTER 2

COMPUTATIONAL DEVELOPMENTS IN POROELASTIC MEDIA

2.1 Introduction

In this thesis, the theory of three-dimensional poroelasticity f i t proposed by Biot (1 94 1, 1955, 1956) is employed to conduct the computational modelling of Iinear, isotropie,

elastic porous media saturated with compressible pore fluids. This Chapter presents the basic partial differential equations goveming the theory of poroelasticity. The essential features of the classical theories of poroelasticity are reviewed with an emphasis on illustrating the physical and mathematical features of Biot's theory of poroeianicity. ï h e basic materiai parameters descnbing the behaviour of proelastic media saturated with compressible fluids are also discussed.

The second part of this Chapter presents a brief review of fuiite element method and boundary element method for the solution of equations goveming poroelasticity. The focus is on illustrating Galerkin's finite element procedure which is used to fornulate the computational scheme employed in this research. The formulations are presented for both plane strain and axisymmetric conditions. The eight-noded plane isoparametric element used in computational modelling will also be discussed.

A Cartesian tensor notation will be used in the presentation with Einstein's summation

convention applied to repeated indices. In addition, the following sign convention is adopted: al1 normal stresses in the solid skeleton and Buid pressures in the pores that are tensile, are regarded as positive. Shear stresses follow the conventional sign convention used in solid mechanics and geomechanics (Fung, 1965;Davis and Selvadurai, 1996). 2.2 Classical Theories of Poroelasticity

The treatment of h s t mathematical formulation and analysis of consolidation behaviour

of soils is due to Temghi (1923). This investigation also introduced the concept of the effective stress principle in the study of geomatenals. Terzaghi (1923) proposed a

hindarnental approach to the study of a Mly saturated soil and developed the onedimensiond theory of soil consolidation by using a mode1 of a porous medium which experiences small deformations. Terzaghi's theory of one-dimensional consolidation has been widely applied to the anaiysis of practical problems and has become a standard procedure in geotechnical analysis. Many conventionai methods for predicting the magnitude and progress of foundation settlements are based on this theory (Hm, 1966; Lambe and Whitman, 1969).

A generalization of Terzaghi's theory to include three-dimensional effects was suggested

by Rendulic (1936). It was based on the assumption that during the process of soil consolidation, the surnmation of total normal stresses in the soil medium remains constant with time. In other words, the first stress invariant I,=a,,+a~+q, is assumed to remain constant with the dissipation of excess pore pressures. This hypothesis has led to

the Terzaglu-Rendulic theory or the pseudo three-dimensional theory of soil consolidation. In this approach, the total stress field in the soil medium is treated independently and usually accomplished by assuming that it is govemed by time independent elastic response characteristics of an elastic medium (see e.g. Zaretskii, 1972).

One of the serious drawbacks of these two theories is the absence of the correct form of coupling between the defonnation of soil skeleton and the flow of pore fluids. As a direct result of this deficiency, some special feahires of the consolidation process, such as the Mandel-Cryer effect (Mandel, 1953; Cryer, 1963) do not appear in these two uncoupled theories. This aspect of the correct form of coupling between the mechanical deformation and fluid flow process is a major development of Biot's theory of three dimensional soi1

consolidation (Schiffinan et al., 1969). The theory of threedimensional linear poroelasticity for a saturated medium was formulated by Biot (1941, 1955, 1956) to mode1 more realistically the mechanical behaviour of sanirated soils or rocks. In this theory, the soil skeleton is modelled as deformable, linear, elastic, porous medium saturated with a compressible fluid. A set of partial differential equations was fomulated by Biot (1941, 1955, 1956) to descnbe the coupled mechanical behaviour of saturated porous media. Biot's theory of poroelasticity

results in a completely self-consistent bounday and initial value problem. Biot's theory inuinsically takes into account the time dependent interaction of the soil skeleton and pore fluids (i.e. the coupling between the deformation of the solid phase and the flow of pore fluids). The coupled mechanical state is described by mechanical variables (effective stresses and excess pore pressure) and kinematics variables (displacements and velocities) within each phase. The mechanicd and kinematics variables are tirnedependent. Four of these variables are independent and represent the coupling characteristics of the goveming equations. They could be regarded as the three displacement components and the excess pore pressure. Consequently, it is necessary and suficient to adopt four independent boundary conditions (three for solid phase and one for pore fluid phase) to set up a proper boundary and initial value problem for a saturated porous medium. These four independent boundary conditions correctly match the

situations encountered with many practical problems associated with saturated soils and rocks.

Biot's theory is sufficiently general in that both the Terzaghi's classical one-dimensional consolidation theory and the Terzagh-Rendulic theory can be recovered as special cases.

2.3 Governing Equations of Poroelasticity The linear constitutive equations governing the quasi-static response of a poroelastic medium, which consists of a porous isotropie elastic soil skeleton saturated nith a

compressible pore fluid take the forms

( 2 .la) (2.1 b)

where a, is the total stress tensor; p is the pore fluid pressure; 5, is the volumetric nrain in the compressible pore fluid; v is Poisson's ratio and p is shear modulus applicable to the porous fabric, 6, is Kronecker's delta function (=1 if i=j; =O if i+j]. Also in (Eq. 2.1 a) E,

is the soil skeleton strain tensor which, for mal1 deformations, is defmed by

where ui are the displacement components, and a comma denotes a partial denvative with respect to the spatial variables. The material parameters a and

P which defme

respectively, the compressibility of the pore fluid and the compressibiiity of the soil fabric are as follows:

where vu is the undrained Poisson's ratio, and B is pore pressure pararneter introduced by

Skempton (1954). The pararneter B is defined as the ratio of the induced pore water

pressure to the variation in total isotropie stress, measured under undrained conditions. This representation of constitutive goveming equations was first proposed by Rice and CIeary (1 976).

In contrat

to

Terzaghi's defuiition of the effective stresses, which is void of

compressibility of the system, in Biot's theory the definition of the effective stress tensor takes into account the compressibility effects. The effective stresses a$ in the porous skeleton are given by

In the absence of body forces and dynamic effects, the quasi-static equations of equilibrium for the entire fluid saturated porous medium takes the fom:

The fluid transport within the pores of the medium is govemed by Darcyoslaw which c m be written as

where ui is the specific discharge vector in the pore fluid and u=Wy, in which k is the coefficient of hydraulic conductivity and y , is the unit weight of pore fluid. The equation

of continuity associated with quasi-static fluid flow is

The bounds for the five constitutive parameten p. v , vu, B, and

K

goveming the

poroelastic behaviour of fluid saturated materials can be obtained by considering requirements for a positive definite strain energy potential (Rice and Clearly, 1976). It c m be shown that the material parameters should satisS the following thermodynamical

constraints:

When dealing with most geomatenals with an elastic response, the lower limit of -1 for v

and vu c m be replaced by zero @esai and Siriwardane, 1984 ; Davis and Selvadurai, 1996). Consequently, we have

In special case of poroelastic material saturated with an incompressible pore fluid, we have

The resulting equations of equilibrium for a poroelastic medium as introduced by Biot (194 1, 1955, 1956) and reformulated in more physically relevant variables by Rice and Cleary (1 976),can be wrinen in terms of the displacements and pore pressure as

The equation (2.9a) corresponds to the displacement equation of equilibrium for linear, isotropie elasticity together with the coupling term of stresses induced by the pore fluid.

The equation (2.9b) is a diffusion-type equation for the pore fluid which includes a

coupling terni to account for deformations of the soi1 skeleton.

The following limiting conditions can be recovered fiom the generalized results: 1. The undrained behaviour îmmediately d e r the loading corresponds to

cu=O,

consequently the equation (2.9a) reduces to the usual elasticity equation with undrained Poisson's ratio vu. i.e.

2. In achieving the &Ily drained state, the excess pore pressure dissipates completely (i.e. p+O) and the equation (2.9a) reduces to usual elasticity equation with drained values of

the elastic constants: p and v . Le.

3. For most naturally occuning geomateriais, the sanirated porous medium can be modelled as a medium consisting of incompressible solid particles and an incompressible pore fluid. In this case

pv2u, +

(1 -

Ci

<, =

EU,

a = 1, and P = m and the goveming equations reduce to:

- p,, = O

Z V ) ~ ~ J

For a well posed problem, boundary conditions and initial conditions on the variables ÿ, p and/or their derivatives can be prescribed.

2.4 Computational Modelling of Poroelastic Media

Numerical methods such as the finite element technique and boundary element technique have been successfblly applied in the solution of problems in poroelasticity which c m

accommodate cornplex geometries and boundary conditions. The application of f ~ t e element procedures to the solution of problems in poroelasticity is given by a number of

researchers Uicluding Sandhu and Wilson (1969) and Ghaboussi and Wilson (1973), who were the first to apply finite element procedures to examine soi1 consolidation problems.

Sirnilarly, the application of boundary element procedures is documented by Cleary (1977) and Cheng and Liggett (1984a).

2.4.1 Finite Element Method

In rhe finite element technique, the domain of interest is subdivided into discrete finite

elements. The elements are comected at nodal points and continuity of displacement and pore pressure fields are enforced at the element boundaries. The values of the field

variables within the elements are interpolated by polynomials of their nodal values. The goveming equations of poroelasticity can then be approximated into a synem of linear matrix equations by application of either the Galerkin technique (Sandhu and Wilson,

1969) or the variational principle (Ghaboussi and Wilson, 1973). The investigation o f

different spatial interpolation schemes and various temporal approximations is given by Booker and SmalI (1975) and Sandhu et al. (1977). Ghaboussi and Wilson (1973) and Booker and Small (1975) have developed finite element procedures for the analysis of problems associated with surface loading of semiinfite media. Selvadurai and Gopal (1986) and Schrefler and Simoni (1987) have used the finite element method to investigate the consolidation behaviour of media of infite

extent, which are modeled by appeal to special infinite elements.

Finite element methods have a greater appeal to engineers and are widely used in engineering applications in dealing particularly with the the-dependent problems and problems involving nonlinear material bebaviour. We adopt a Galerkin f ~ t element e procedure to formulate the computational procedures for the solution of the various types of poroelasticity problems discussed in Chapter 1. 2.4.2 Boundary Elemeat Method

The boundary integral equation method is a very powerhl technique for the solution of boundary value problems for unknown displacements, pore pressures. and tractions on the boundary of domain. This approach can also provide solutions for intemal field variables, although die finite element analysis is considered to be much more efficient for this purpose.

The boundary element method for non-dissipative systems commences with the use of Betti's reciprocal theorem (see e.g. Davis and Selvadurai, 1996) which relates the work done by two different loadings on the same body. The governing equations of poroelasticity are formulated in a Laplace transform space and then integrated over the boundary of the domain. Betti's reciprocai theorem is then applied

to

the boundary

conditions to write the goveming integral equations in terms of boundary unknoms (in contrast to finite element methods). The boundary is discretized into elements using polynomial approximation of the boundary geometry, displacements and tractions. The resulting integral equations are nurnerically integrated to give a set of linear matrix equations. Finally, the Laplace transform is nurnerically inverted in order to obtain the solution in physical space. Cleary (1977), Banej e e and Butterfield (1981), Kuroki and Onishi (1 982), Cheng and

Liggett (1 984a), and Domingua (1992) have developed boundary element procedures to problems associated with soi1 consolidation.

2.5 Finite Element Formulations The standad Galerkin (1915) finite element procedure (Zienkiewicz, 1979) can be applied to approximate the equations governing Biot's theory of poroelasticity (Eqs. 2.9q b) and reduce them to a system of linear matrix equations. The details of these procedures

are well documented by Sandhu and Wilson (1969), Aboutit et al. (1985)' Schrefler and Sirnoni (1 987), Lewis and Schrefler (1987), Selvadurai and Kapurapu (1989) and more recently by Selvadurai and Nguyen (1993) in comection with the îïnite element modelling of thermal consolidation of sparsely jointed porous media. A bief description of the Gaierkin technique and its application to poroeiasticity problems is given in the following sections.

2.5.1 Galerkin Weighted Residual Method The Galerich method has been widely applied to the fuite elernent formulation of consolidation problems (Lewis and Schrefier 1987). The Gaierkin technique is a speciai case of the general method of weighted residuals. The method provides an approximate solution to the system of partial differential equations of following form:

where L is an operator; 4=$(x} is the unknown solution (function of the coordinates x. ; ,) f=f(x) is a known fùnction of x.; and R is the bounded domain where the equatioa is J

defined. Consider n discrete points (nodes) of the domain R where @=$ An approximate value of

equation:

at each point I.

4 at any arbitrary point of R is assumed by the folloxing

where N,=.V(x) are interpolation or shape fùnctions (usually polynomials). Since +Û(x) 1 J

is not the tnie solution, there is an error (or a residuai) Re associated with the approximation procedure which can be obtained by replacing the value of

4

by û in

Equation 2.13 as follows:

The method of weighted residuals minirnizes the error between the m e solution $I and the approxirnate solution U by reducing the weighted average value of residual Re over the

domain R to zero. The nodal values $, are detemiined by solving the following constraint equations:

where fi; are weighting fùnctions which are functions of the coordinates x,. In the J

Galerkin method the weighting fùnctions W, are assurned to be identical to shape huictions .V,, When al1 nodal values 4, are obtained, Equation 2.14 can be used to obtain the approsimate solution for any arbitrary point in domain R.

2.5.2 Finite Element Approximations

In the solution of an initial boundary value problem, the goveming equations 2.9a and

2.9b and appropriate initial conditions need to be satisfied within the domain R and appropriate boundary conditions should be satisfied on the boundary B of the domain.

Using Gaierkin's technique, the goveming equations can be transformed into matrix

equations where the unknowns are the nodal displacements and pore pressures. Letting u, (i=l, 2, 3; b1,...,N) be the nodal displacements for N discrete points in R, and pK (K=l,

..., n) the nodal pore pressures in n nodes at an arbitrary t h e t. The displacement vector

and pore pressure for any arbitrary point with coordinates x.J in the domain R are approximated by the following relations:

where

Ni and N [ correspond to shape functions for displacernent and pore pressure

fields; J=1, ..., !V and K=l,

..., n; u, is the displacement of the soi1 skeleton at node J in

the ifh direction. The indices in capital letiers (e.g. Jand K ) refer to nodal values, whereas the indices in small leîters (e.g. i andfi refer to coordinate directions. Also the surnrnation

convention is adopted. In general N; and N,P cm be different but both N I and N: must exhibit CO continuity.

The application of a Galerkin procedure to the governing equations (see e.g. Sandhu and Wilson, 1969; Lewis and Schrefler, 1987; and Selvadurai and Nguyen, 1995) gives rise to the following discretized forms of the equations goveming poroelastic media:

K

where

K = stiffiess rnatrix of the soi1 skeieton; C = stifniess matrix due to interaction between soi1 and the pore fluid;

E = compressibility matrix of fluid; H = permeability ma&; F = force vectors due to extemal tractions, body forces and flow field; u,, p, = nodal displacements and pore pressures at tirne i;

At = time increment and y is a time integration constant;

and (

denotes the transpose.

When ut and p, are known, the solution of Equation 2.19 resulu in ut,, and pl+,,. Expression for the matrices K, C and etc. are given in the Appendix A. The tirne integration parameter y varies between O and 1. The criteria for the stability of the integration scheme given by Booker and Small (1975) require that 7 2112. According to Lewis and Schrefler (1987) and Selvadurai and Nguyen (1995), the stability of solution can generally be achieved by selecting values of y close to unity.

Since for many diffusion-type phenomena, the time variation of displacement and pore pressure field variables in pore fiuid diffusion process is exponential in nature, Sandhu and Wilson (1969) have used a logarithmic interpolation function in the tirne solution of

Equation 2.19. This requires the assemblage and triangularization processes of the general matrices at each tirne step Aî of the analysis. In order to rninimize the computational effort and irnprove the numerical efficiency, the t h e sep Af can be logarithmically altered after a selected number of time steps. 2.5.3 lnstantaneous Poroelastic Respoase

The undrained elastic behaviour of saturated materials resulting imrnediately afier the application of extemal loads (when there is not enough time for pore fluid migration within the porous medium) can be incompressible. However, the soi1 particles and pore fluid of most porous sahuated geomaterials are slightly compressible. For this reason, it is usehl to retain the constitutive modelling such that the undrained behaviour is compressible.

The numerical trament of incompressible elastic matends was f k t examined by

Hemnann (1965) using a displacement constraint formulation. Christian (1 968) modified

Figure 2.1 Solid isoparametric element.

the Hemnann formulation to examine the undrained behaviour of consolidating soils. Naylor (1974) has s h o w that the conventional fhte element method c m be applied to nearly incompressible materials ( ~ 0 . 4 9 9 )when a reduced (2x2) integration scheme is used in the analysis.

In this thesis, the conventional finite element method is used to examine the instantaneous behaviour of poroelastic geomaterials which c m exhibit initially compressible elastic behaviour. 2.5.4 Finite Element Discretizations

Following the standard finite element procedures, the domain R is discretized into certain

number of subdornains called finite elements. The element chosen to represent the intact region of the poroelastic medium is the eight-noded plane isoparametric element where the displacements within the element are interpolated as functions of the 8 nodes, whereas

the pore pressures are interpolated as a fùnction of only the four corner nodes i, k, m, and O

(Figure 2.1). This type of element has k e n used to mode1 the consolidation behaviour

of poroelastic media by several researchers including Aboutit et al. (1 985). Smith and

Figure 2.2 Tirne variation of prescribed displacement and traction boundary conditions. Griffiths (1 988), Selvadurai and Karpurapu (1 989), and Selvadurai and Nguyen (1 995).

Aboutit et al. (1985) have s h o w that this type of element is less prone to spatial oscillations in the solution obtained for the pore pressure as opposed to the 8-noded element type where al1 degrees of fkeedorn are calculated at al1 the nodes. The physical

explanation of this spatial oscillations can be argued as follows: fluid pressure has the sarne dimension as the stress. Strain is directly related to stress via the elastic constitutive parameten (e.g. p and v). Since serain is expressed in tems of spatial derivatives of

displacernents, the polynomials used as interpolation functions for fluid pressure should, for consinency. be one order Iower than the ones used for displacernents. 2.5.5 Time Variation of Bounda y Conditions

To examine the effect of time variation and rate of Ioading on the poroeiastic behaviour of geomaterials, the variation of prescribed displacement and traction boundary conditions with time can be taken into account in the finite element formulations. In this

study, a time variation in fonn of a ramp function (Figure 2.2) for the imposed displacements and tractions is considered as follows:

The time-dependent response of poroelastic media associated with imposed boundary conditions given by Equation 2.20 corresponds to the settlement consolidation under a prescribed load, or load relaxation under a prescribed displacement. The limiting value of

the rate of loading parameter, when m=- (Figure 2.2), corresponds to a Heaviside step function (i.e.At)=H(t))where the entire load is applied instantaneously.

VERIFICATION OF THE FINITE ELEMENT CODE

Venfication of a cornputer code is the process whereby the reliability and accuracy of the numerical solutions to the fundamental equations are docurnented. In order to ver@ the code, the results derived frorn the computationai mode1 are compared with a series of

benchmark problems, which are primarily analytical solutions. The analytical solutions

for one-dimensional consolidation are given by Temghi (1923), and analytical results for consolidation of a poroelastic sphere are given by Gibson et al. (1963). The anaiytical solutions for two-dimensional plane strain or axisyrnrnetric consolidation of half space and layer regions are also given by McNamee and Gibson (1960b), Gibson et al. (1970),

Chiarella and Booker (1975), Selvadurai and Yue (1994) and Yue and Selvadurai (1995). In this Chapter we present a venfication of the finite element procedure developed in comection with the research, by cornparison with certain anaiytical solutions.

3.1 One-Dimensional ConsoIidation

Terzaghi (1923) presented an analytical solution for the problem of the one-dimensional consolidation of a soil column of depth H which rests on an impermeable base. The medium is subjected to a total stress

with time variation in the f o m of a Heaviside

step h c t i o n . The finite element mesh and the boundary conditions for this problem are

s h o w in Figure 3.1.

The classical anaiflicai solution for the evolution of pore pressure p within the soil column which satisfies the drainage and displacement boundary conditions indicated in Figure 3.1, takes the fom:

where M = ( x / 2 ) ( 2 i+ l),

i = 1,2,. .. and

Similarly, the analytical solution for time variation of the surface seulement ( w ) takes the form:

To conduct the comparisons, the following matenal properties are useci in the finite

element andysis:

Figure 3.1 Finite element discretization of one-dimensional consolidation.

The compressive total stress applied at the surface is assumed to be a*= 10x10~kPa; and the depth H of soi1 deposit is assumed to be 100 m. Figw 3.2 shows that the finite

element results for the the-dependent variation of pore pressure agrees well with the analytical solution. Figure 3.3 illustrates the numerical and analytical results of normalized surface senlement which are in very good agreement (maximum discrepancy

is 5%). The above numerical results were obtained with a value of the integration constant y=0.875. For lower values of y, some instability occun in computed pore pressure at small values of the nondimensional time factor T(R0.0001).

Numencal - r/H--0.597 a

Numerical - r/H=O.1 15

- Andytical - rlH=0.597 - - Analytical - dH=O. 1 15

0.0001

0.001

o. 1

0.01

I

10

T

Figure 3.2 Pore pressure evolution at different depths.

Figure 3.3 Time variation of d a c e senlement.

3.2 Consoüdation of a Poroelastic Sphere

In the process of soi1 consolidation, some regions of materiai experience a delay in the development of a peak in the time-dependent response of the pore pressure. The above phenornenon is called the Mandel-Cryer effect and was mathematically demonstrated by Mandel (1950, 1953) for a cubic body uniaxially loaded under plane strain conditions, and Cryer (1963) for a sphere subjected to a uniforni pressure at its fiee drainage surface.

The Mandel-Cryer effect has also been experirnentally observed by Gibson et al. (1963)

and Verruijt (1963) in saturated geomaterials such as clay.

For the example of a

uniformly loaded sphere, the Mandel-Cryer effect can be physically explained as follows (Cryer, 1963). At early times in the consolidation process, almost al1 of the volume change occurs near the surface where fluid draining occurs. The region near the surface will tend to contract resulting in squeezing of the central regions. Consequently, the total radial stress in the central regions will increase. As there is little volume change in the

central regions, the effective stress remains constant and therefore the pore pressure will rise. The consolidation problem of a poroelastic sphere of radius a, subjected to a uniform

total stress a. at its &e draining sunace is considered. The analyticai solutions for the time-dependent variation of pore pressures to this problem are given by Gibson et al. (1963). Figure 3.4 shows the finite element discretization of the sphere and associated boundary conditions. Figure 3.5 illustrates the analytical and numerical results of the pore pressure variation at the center of a sphere for difTerent values of Poisson's ratio v. The general agreement is considered to be satisfactory.

3.3 Plane Stnin and Axisymmetric Consolidation McNamee and Gibson (1960b)developed analytical solutions for the plane strain and the axisymmeûic problems related to consolidation of a poroelastic half-space subjected to a unifom total normal stress oo.The loaded region for plane strain conditions is a strip

a

impermeable

a

--

Figure 3.4 Finite element discretization of poroelastic sphere.

Figure 3.5 Pore pressure evolution at the center of poroelastic sphere.

load of width 2a, and for axisymmetric conditions the loaded region is a flexible circular region of diameter 20 (Figure 3.6). McNamee and Gibson (1960b)developed integral

expressions for the displacements, stresses, and pore pressures. These i n f i t e intepis can be numerically evaluated by ushg a suitable quadrature scheme. For the special case of Poisson's ratio v=O, the expressions c m be simplified to:

where

-

X=xla (rla for avisymmetric case), Z= zla and

K ( X , 5 ) = 1 n;

COS(-)

sin 6 for plane strain

erf is the error function defined by:

JO

and J I correspond to Bessel fùnctions of the first kind of order O and 1, respectively.

The finite element discretization and the boundary conditions for this problem are shown

in Figure 3.6. The following parameters are used in analysis:

The numerical and analytical results for pore pressure evolution are compared in Figures 3.7 for plane strain and axisymmetnc conditions. In the numerical computations, a value of y=0.875 was used. The agreement is quite good considering the influence of finite

domain used in finite element modelling of a poroelastic half-space. The f ~ t element e method tends to overpredict the values of pore pressure with the boundary conditions sh0u-n in Figure 3.6. Typically, at the time corresponding to the peak of the time-

dependent pore pressure, the numencal value for the pore pressure is 3% higher than the analytical value; at later times (T=IO),this overprediction in numericd result reaches a value of 10%. A second analysis, with a zero pore pressure condition specified at the

right hand side boundary, was performed. That resulted in an underprediction of the

results obtained by finite element method, with typical absolute values similar to the former case. For both plane strain and axisymmetric conditions, the Mandel-Cryer effect is manifest at al1 points k i n g considered. The integrals occming in Equation (3.4) are numerically evaluated using the symbolic manipulation code MATHEMATKA (Wolfram Research

Inc., 1993, V.2.2).

Gibson et ut. (1970) examined the consolidation problem of a poroelastic layer laying on a rigid base. They developed integrai expressions for time variation of the surface

senlement and presented the results for the tirne-dependent variation of the degree of settlement consolidation U(), defined as follows

Degree of consolidation provides a measme of the extent to which the time-dependent part of consoiichtion has been attained. Figure 3.8 illustrates cornparison of the numerical

Free drainage @=O)

impermeable

b

impermeable Point A:

A

,-a

Point B:x=O

~ 2 . 7 5 ~

Point C:x=a

~=2.75a

Figure 3.6 Finite element discretization of the poroelastic hdf-space under plane strain or axisymmetric consolidation. and analyticai results for the degree of consolidation (h/u=lO) in plane strain and axisymmetric conditions. The maxirnum discrepancy is 10Y0. The agreement is

considered to be satisfactory for purposes of geotechnical applications.

P.A (Numerical)

- P. B (Analytical) O

a-

-...O

P. B (Numerical)

P. C (Nurnerical) P. C (Analytical)

a) Plane strain

-

P. A (Numerical) P. B (Analytical)

0.4 O

0.3

P -

P. B (Numerical) P. C (Numerical)

0 0

0.2

o. I O

Figure 3.7 Tirnedependent pore pressure evolution.

Numencal

- Analytical

a) Plane strain

Numencal

- Analytical

b) kvisymmetric Figure 3.8 Degree of consolidation for proelastic layer.

FRACTURE MECEANICS OF POROELASTIC MEDIA

4.1 Introduction

.4pplications of the theory of poroelasticity have largely focused on problems related to conventional geomechanical applications where the emphasis has been on the evaluation of deformations and consolidation rates in such media subjected primarily to extemal loads. As the applications of the theory diversi@, attention needs to be focused on other aspects of importance. With prous brinle geomaterials which are saturated with fluids,

the initiation and extension of cracks are recognized as

ui

area of both practical and

fundamentai interest. This Chapter focuses on the finite element rnodelling of the fracture phenomena in pon>elmïc media The class of transient stationary cracks, where the migration of pore fluid into the m c k tip governs the time-dependent behaviour of stable

cracks with constant iength, is examhed for both plane straùi and axisymmetric

situations.

In the context of analytical studies the treatment of cracks occurring in poroelastic geomaterials has received ody l u t e d attention, The problem of the shear fracture in a poroelastic material was investigated by Rudnicki (1985) in connection with fault initiation and propagation in porous rock. This problem was examined in view of its potential use in the study of earthquake mechanisms at the source location. The class of

transient fracture problems, where there is the influence of pore fluid processes, are important in connection with enhanced energy resources recovery by the process of hydraulic fracturing (Detomay and Cheng, 1991). In hydraulic fkacturing techniques,

during the initiation of hcture by the pressurization of the fluid, the pore fluid tends to leak out into the porous formation and induce a volume expansion of material around the crack tip region leading to partial closure of the fracture. Field evidence (see e.g. Smith, 1985; and Nierode, 1985) suggests that the injection pressure required to initiate the

fracture in poroelastic media is higher than those predicted by classical eiasticity models which neglect the influence of poroelasticity. Recently Craster and Atkinson (199 1, 1996) have examined the problem of stationary semi-infinite cracks embedded in poroelastic

media These studies take into account the tirne-dependent behaviour of the crack in the

presence of variable pore pressure boundary conditions at the faces of the crack and assess their potential influence on the pore pressure field and the effective stress intensity

factors at the crack tip. The numencal modelling of stationary hctures and hcnire extension in geo1ogica.i media can be approached by adopting either finite element techniques or boundary

element techniques. Accounts of these developments together with extensive references to literahire in these areas are given by Zienkiewicz (1977), Brebbia et al. (1984), Brebbia and Aliabadi (1 !IN),and Atluri (1991). The computational modelling of transient hcture

problems in poroelastic media has been largeiy examined in relation to the hydraulic

fractimng techniques. Boone and Ingraffea (1990) and Iograffea and Boone (1988) have examined the two-dimensional propagation of a plane crack driven by the pore fluid

pressures in a poroelastic matenal. They developed a crack extension mode1 based on the separation of the interface elements placed dong an apriori kwwn crack path when the

maximum tensile stresses at the crack front reach a cnticd value. This procedure ignores the kacture resistance of material which c m be significant in fracture phenornena of poroelastic media in the absence of relatively high in situ stresses. Vandamme et al. (1 989) have examined, via a displacement discontinuity method, the transient plane strain behaviour of a stationary poroelastic fracture subjected to a sudden constant pore fluid pressure. They decomposed the fluid pressure loading into prescribed pore pressure

boundary conditions and normal traction boundary conditions on the crack faces and examined the time-dependent behaviour of the crack separately for each type of boundary condition. Detournay and Cheng (1991) have examined similar problems in more detail to evaluate the design parameters, such as the time-dependent variation of the crack opening stress intensity factor and the pore fluid volume leak-out into the porous

medium, related to the process of hydraulic frachiring by using a similar modelling technique. These computational models can be extended to cover the geomechanics problems where the crack faces cm be subjected to variable traction, displacement and pore pressure boundary conditions in poroelastic media. The numerical procedure developed in this Chapter considers the interaction eflects between the migration of pore

fluid and the deformation of a porous fabnc on the transient behavior of poroelastic cracks which are subjected io general loading modes.

In this thesis, a finite element procedure is developed to examine b o t . axisymmetric and two-dirnensional plane straui problems of fracture mechanics associated with poroelastic media. The computationai mode1 takes into account the stress singularity of the effective stresses at the crack tip. The computationai results are verified by cornparison with the

analytical solutions given by Craster and Atkinson (1991, 1996) for the cases of permeable and impermeable semi-ulfuiite cracks embedded in poroelastic media of

infinite extent.

The numericd procedure is also utilized to examine the axisymmeûic surface indentation of a poroeiastic half-spacc region which is weakened either by a cylindrical crack or a penny-shaped crack. The results illustrate the influence of the extent of hcture and the

pore pressure boundary conditions on the various surfaces. on the time-dependent evolution of the stress intensity factors and the tirnedependent consolidation settlement of the axisymmetnc smooth @id indentor. 4.2 Fracture Mechanics Concepts

The linear elastic fracture mechanics (LEFM) is based on the theones of fracture onginally proposed by Griffith (1921, 1924). The classical rheory of elastic hachue mechanics is based on the linear theory of elasticity. This theop deals with sharp cracks and assumes that al1 the nonelastic phenomena associated mith the hcture process takes

place only at the crack tip. This is a suitable assumption when dealing with cracks which

are small in cornparison to dimensions of elastic body, which is the case in most problems in geomechanics. In this section we will present a brief review of some of basic concepts and definitions employed in linear hcture mechanics.

The theory of linear elastic fracture mechanics originated fiom classical stress analysis of the elliptical fiaw problem examined by Inglis (1913). G S t b (1921, 1924) then

extended this snidy to examine the crack problem which defmes the conditions necessary

for the extension of the crack. He developed a fracture theop based on the concept of energy balance and proposed that a flaw becomes unstable when the elastic strain energy

change, which results from an increment of crack growlh, is larger than the surface energy (i.e. resistance) of matenal. The assumption of linear elastic behaviour is central

to the developments proposed by Griffith (1921, 1924). The concepts proposed by

Griffith were extended by Irwin (1957, 1958) to examine the behaviour of the local effects at the tip of sharp cracks. Irwin (1 957, 1958) adoptai the analytical technique developed by Westergaard (1938) in analysis of stresses and displacements ahead of a

s h q crack in elastic materials, and showed that the crack tip solutions could be described by a single constant. This parameter which characmixs the displacement and stress fields

near the crack tip of hctures, later became homn as the stress intemity

fictor. Williams (1957) applied a different technique to derive the d y t i d solutions for

in-plane opening

mode 1

in-plane shearing mode U

anti-plane shearing mode i l

Figure 4.1 Fundamental modes of deformation at the crack tip region. the stress state at the crack tip which were essentially identicai to those obtained by Irwin (1 957, t 958).

The assumption of h e a r elasticity in elastic fisictue mechanics allow-s us to identify various modes of deformation of the crack tip region which can initiate the growth of cracks in proelastic materials under the influence of effective stress fields. Invin (1957, 1958) observed that there are in generaf three independent modes of deformation at the

crack tip region for elastic materials. These three modes of hcture can in general be

classified as mode I (in-plane opening of crack), mode II (in-plane sheariog of crack) and mode III (mti-plane shearhg of crack). Figure 4.1 illustrates these fundamental modes of deformation at the cnck tip region. The influence of al1 extemal stress nates can be represented as a linear combination of these three modes of deformation at the crack tip. Only the crack opening mode 1 and crack shearing mode II need to be considered in fhcture analysis of plane elasticity problems. Invin (1957, 1958) showed that the near

crack tip asymptotic field of stresses a, and displacements 4 in an isotropie elastic

material under plane sûain conditiom take the followîng foms:

r and 0 are polar coordinate reference centered upon the crack tip (Figure 4.1). 1;,, g,, hi,

and w, are functions of 0 and depend only on the mode of crack deformation. The indices

KI,are referred to as the stress intensity factors associated with the two in-plane modes of fracture. The mess intensity i and j refer to x and y directions. The consta~tsKIand

factors depend on the Ioading conditions and characterize the magnitude of the crack tip stress field. If the crack lies on the axis x (Figure 4 4 , the following functions can be

denved for the crack opening mode I deformation at the crack tip in plane main problems

(see e.g. Broek, 1982)

0

38 2

(1 T sin - sin -)cos(?)

2

0

The stress components related to the plane ofû=û at the crack tip region associated with the two modes of deformation, take the following forms

Equation 4.3 illustrates that, for elastic materials, the stress field possesses a singularity of order i l R at the crack tip. It is assumed in these formulations that the angle between

the faces of crack is zero. The stress field at the crack tip shows a general singularity of order rL1for large crack angles.

For poroelastic materiai, the behaviour of stationary cracks is govemed by the time-

dependent pore pressure diffusion process at the crack tip region which depends on the pore pressure boundary conditions at the crack faces. The main energy due to application

of extemal tractions in poroelastic media is reduced by the migration of pore fluid into the crack tip. As the pore pressure drainage takes place near the crack tip region. the

stress intensity factors Vary from the limiting case of elastic undrained behaviour (with Poisson's ratio vJ to the fùlly drained bebaviour (with Poisson's ratio v) in poroelastic materials. This tirnedependent behaviour of cracks is a characteristic feanire of fractures in poroelastic media. Through a boundary Iayer analysis of poroelastic crack problems, Sirnom (1977) has shown that the stress singularity of order

is preserved for the effective stress field at

the crack tip in poroelastic media. 'This stress singdarity is maintained throughout the tirne-dependent poroelastic bebaviour of saturated materials. Craster and Atkinson ( 1 99 1) have shown that the pore pressure behaviour at the crack tip as r+ O is not spatially

singular for transient poroelastic hcture problems. They have show that the pore pressure gradients at the crack tip are, however, singular for the permeable crack problems. The following argument can be used to show that the pore pressure exhibits temporal singularity at the crack tip. The pore pressure can be related

to

the value of

mean total stress (ad3) through Skempton's pore pressure coefficient B in undrained behaviour (H')

of poroelastic materials. Since the stress field at the crack tip shows

singular behavior, the pore pressure field also exhibits singularity in the undrained behaviour of proelastic matenal correspondhg to time Hl+.However, as the pore fluid drainage takes place at the crack tip region, the pore pressure field becomes reguiar.

4.2.1 ModeMing of Crack Tip Behaviour

In b i t e element modelliog of continuous media, tbere are always computational errors associated with the f ~ t element e discretization of domain of interest. The approximate solutions of the finite element method can converge to the exact solution obtained by the

classical elasticity. The satisfaction of displacement continuity at inter-element boudaries and the ability to mode1 constant strain and rigid-body displacements ensures that the finite element solutions converge to the tnie solution monotonically with

progressive mesh refïnement, as the size of finite elements approach to zero (Zienkievjicz, 1977). An indication of the discretization enor associated with finite element solution can be obtained by the following

where uoand

u' correspond to the exact solution and the finite element approximation of

problern in domain R. Tong and Pian (1973) have show that the rate of convergence in the f ~ t element e method for a plane elasticity problem with a stress singularity of rLi

can be established by following

e(x) 5 ch'

(4.5)

where h is the maximum size of the elements and c is a positive constant which depends on the elastic constants, geomeûy of mesh, and behaviour of u near the point of singularity. This indicates that the rate of convergence in the presence of a singularity is controlled only by the nature of the singularity and is independent of order of

polynomiais used in the interpolation functions. If there is no stress singularity, the convergence rate in displacements will be of order h2 for quadratic isopararnetnc elements (Tong and Pian, 1973). The presence of a stress singulanty slows down the convergence rate and may result in loss of accuracy when the piecewise polynomial interpolation f'wictionsare employed

b

-

a

7

6

A

5$

~,,~~*

crack tip 4

,2

1.1

crack tip -

-u 4 -

3

4

Figure 4.2 Quarter point singuia.elements. in fmite element aoalysis. in order to improve the rate of convergence, the interpolation

fùnctions should include tenns that c m account for the proper singularity. Tong and Pian ( 1 973) suggested that the special crack tip elements should be employed to minimize the

discretization enor given by Equation 4.4. Henshell and Shaw (1975) and Barsourn (1976) independently proposed the use of crack tip singular elements for the f ~ t element e modelling of linear elastic media which contain defects such as cracks. The procedures adopted in these analyses were to shift the mid-side nodes of an eight aoded isoparametric element to their quarter points. This

relocation of mid-side nodes gives rise to the required singularity in the m;iin field and consequendy a

116 stress singularity at the crack tip. For example refemng to Figure

4.2, the basic bctions for the displacement fields associated with the regular isoparametric element dong side 1-3 c m be written in the fom

When the mid-side nodes are shifted to the quarter point locations (i.e. x,=û, x2=t,/4, x3=tl), the expression for x coordinate takes the fom:

and solving for local 5 coordinate

Substihiting Equation 4.8 into Equation 4.6, the displacement variation dong side 1-3 can be written in the fom:

The correspondhg strain in the x-direction takes the fom:

This shows that the quarter point singular element has the required strain singularity of order Y '

along the two edges containhg the crack tip. However, it cm be shown by

using a similar approach that along any ray within the element emanating fiom the crack tip, the

variation is not of the form

fin. This condition can be obtained fiom a

triangular element by collapsing one side of quadratic element onto the crack tip node (Figure 4.2). Again the mid-side nodes close to the crack tip are moved to their quarter points. The latter form of singular crack tip element has s h o w better performance with more accurate results (Barsoum, 1976). The singular crack tip elernents have k e n proved

to satisfy the inter-element compatibility and continuity conditions (Barsoum, 1976). Tfie quarter point crack tip elernent offers its preference over any other special crack tip

elements (such as hybrid elements and etc., Owen and Fawkes, 1983) due to its simplicity and its ready availability in any finite element code. This element has k e n successfùlly utilized to mode1 crack problems in classical elasticity via finite element methods (Ingraffea, 1977% b; Owen and Fawkes, 1983; Murti and Valliappan, 1986) and boundary rlement methods (Cruse and Wilson, 1977; Blandford

et

al., 1981; Smith and Mason,

1982; Selvadurai and Au, 1989; Selvadurai and ten Busschen, 1995).

In this study the crack tip element developed by Henshell and Shaw (1975) and Barsourn (1976) is adopted to mode1 the stress singularity in the effective stresses at the crack tip applicable to the porous skeleton. The pore pressure field near the crack tip is, however, rnodeled by regular quadrilateral isoparametric elements. 4.2.2 Evaluatioa of Stress Intensity Factors

The stress intensity factors at the crack tip can be evaiuated by using the results for the displacement and stress fields in the vicinity of crack tip given in Equations 4.1. It can be shown (see e.g. Shih et d,1976) that the coefficient of the first terni of the analytical solutions for the near tip displacement fields given by Williams (1957) can be related to

the displacements computed via the finite element method at the equivaient location. This

Figure 4.3 Node arrangement for computation of the swss intensity factors.

procedure is referred to as the displacement correlation method. In particular this method utilizes the nodal displacements at four locations A, B, E, D of the crack faces as show in Figure 4.3. The flaw opening mode 1 and flaw shearing mode II stress intensity factors take the fonns

where k, = (3-4v) and L, is the length of the crack-tip element The accuracy of this approach has k e n assessed by Ingraffea (1977a) and Murti and

Valliappan (1 986) in computational hcture mechanics of elastic media. This procedure

is found to give excellent cornparison with known analytical solutions when the following simple guidelines are followed These include: (1) a reduced numerical

integration scheme (2x2) for the crack tip elements; (ii) an aspect ratio of element close to unity; (üi) the size of crack tip element (4,) in the range of 1525% of the crack length;

and (iv) the angle of the singular element at the crack tip Iess than n/2. 4.3 Stationary Poroelastic Fracture

h this section, the transient behaviour of a semi-infite stationary Fracture embedded in a poroelastic medium is examined. Specifically we consider a crack configuration as shown

in Figure 4.4. The problem can be regarded as an extension of a Griffith-Williams andysis of a semi-infinite fracture to include poroelastic effects. The tensile normal traction hundary conditions on fracture surfaces result in a crack opening mode I deformation at the crack tip region. The pore pressure boundary conditions on the faces of the crack c m either be permeable or impermeable. The symmetry of the problem about

the plane F O allows us to consider a single half-space region occupying z >O. The

imperneable pore pressure boundary conditions ahead of the crack tip dong the crack axis (x>O) is considered due to the symrnetry in mode I deformation. The interaction

between the deformation of the porous skeleton and the flow of pore nuid results in a time-dependent behaviour of fractures in poroelastic media. The stress intensity factors which characterize the singular crack tip stress fields are time- dependent, and depend

upon the pore pressure boundary conditions on the ûacture planes. The initial and fmal values of stress intensity factors, however, correspond to the undrained and drained eiastic behaviour of the material.

It can be shown that the time dependent-variation of the stress intensity factors can be computed quite conveniently using the h i t e element technique ernploying the quarter

point singuiar elements.

crack tip

7

Figure 4.4 Finite element discretkation of a semi-infinite crack

in an infinite proelastic medium. 4.3. f Verification Exercises

The numerical procedure for evaluation of the stress intensity factors has k e n verified by appeal to the analytical solutions. We consider the problem of a semi-infinite crack in an isotropie poroelastic medium, the faces of which are subjected to an intemal normal total

stress a,which decays exponentialiy with respect to the spatial variable x. i.e.

where H(r) is the Heaviside s e p fùnction, and a is a constant. Figure 4.4 illustrates the geometry and the fiaite element discretization of the problem and associated bouodary conditions used in the computations. In Figure 4.5 the results derived fiom the Gnite element technique for the t h e variation of the flaw opening a e s s intensity factor K,are compared with the analytical results given by Craster and Atkinson (1996). The

normalized time factor T used in results is given by

where C is the genedized consolidation coefficient for the compressible porous medium defined by

The material parameters v

=

0.3,vu = 0.4 are assumed in analyses and two cases of

permeable and impermeable pore pressure boundary conditions for the crack faces are

considered. The results indicate tbat there is good agreement between analytical results and cornputational results particldarly with respect to the time variation of the stress intensity factor. The maximum discrepancy of 1.5% in the results, which is obtained for

the permeable crack, could be amibuted to the relatively coarse finite element mesh

discretization at the crack tip.

The time variation of the stress intemity factor is due to the migration of pore fluid

through the poroelastic matenal. At large times; as T + a,the region around the crack tip

will be drained and the bebaviour of the poroelastic medium is elastic. At small times the

-

U

5

T (Tirne factor) (a) Permeable crack

#' 0.98

0

z O

Numerical

"

CI

.

œ

- = g g

i-)

T ( T i c factor)

(b) impermeable crack

Figure 4.5 Time dependent behaviour of the flaw opening mode stress intensity factor KI.

movement of fluid towards the crack tip lowen the stress intensity factor. This results in a reduction of available energy to h c t u r e the material. Alkinson and Craster (1991) showed that the pore pressure at the crack tip is singular with time. This can be attributed to the behaviour of the crack tip as a pore pressure "sink".

4.4 Indentation of a Cracked Poroelastic Haif-Space

The earliest application of the theory of poroelasticity to initiate bondary value problems in geornechanics commenced with the classical studies by McNarnee and Gibson (1960a, b) who applied integral transform techniques to the solution of surface loading problems

dealing with a poroelastic half-space region. Further application of anaiytical studies

include the works of Cryer (1963) and Schiffman and Fungaroli (1965). These works were followed by some important developments involving interaction of rigid indentors

and half-space regions (Agbezuge and Deresiewicz, 1974, 1975; Chiarella and Booker,

1975; Szefer and Gaszynski, 1975; Booker and Small, 1975). More recent analytical

research in the area of poroelastic contact and inclusion problerns is due to Selvadurai and

Yue (1994) and Yue and Selvadurai (1994, 1995). This section deals with the application of computational modelling to the study of the

classical indentation problem associated with a fluid sanirated half-space region. The indentation of the surface of an elastic half-space by a smooth flat indentor was first examined by Boussinesq (1885) by applying results of potential theory. The problem was re-examined by Harding and Sneddon (1945) who reduced the axisymmeeic problem to the solution of a pair of dual integral equations. The Iiteratw in contact mechanics contains extensive accounts of the subject of indentation problems and such studies are given by Galin ( 196 1), Lur'e (1W), Selvadurai (1 979), Gladwell (1WO), and Johnson (1 985). In the classical treatment of the axisymmetric contact problem, it is invariably

assumed that the elastic medium is fke of defects, most notably cracks, which are present

in the vicinity of the indentation zone. The primary reason is tbat the presence of defects even in their axisymmetric forms entail a great deai of analytid effort in the solution of

the final foms of the integral equations (Selvadurai, 1997b). The applications of the classical theory of poroelasticity to problems of indentation in the vicinity of a defect

makes the problern unwieldy particuiarly due to the convolution nature of the fmal integral equations governing the poroelasticity problem (Lm and Selvadurai, 1996). For

this reason, recourse must be made to the use of numerical schemes when dealing with problems associated with interaction of indentors and cracks which can develop and extend in the vicinity of the indentation zone.

The analysis focuses on the axisymmetric indentation of a saturated proelastic halfspace which is weakened by a stable crack. The geometry of the crack is considered to be

axisymmeiric and two configurations of the crack are considered; the tint involves a cylindricd crack which is located at a specified radius fiom the ngid indentor (Figure 4.6a), and the second involves a tlat penny shaped CO-axialcrack which is located at a

finite depth hom the rigid indentor (Figure 4.6b). The computational modelling primarily

focuses on the examination of two effects; the first relates to the influence of stress singularity on both thedependent displacement of the indentor and pore fluid pressure at the crack tip, and the second relates to the influence of the defect and pore pressure

boundary conditions on the time-dependent displacements of the indentor. The relative geometry between the dimension of the indentor and the dimension of the crack is taken as a variable. 4.4.1 Cyiindrical Crack

The problem of a deep deposit of a saturated poroelastic geomaterial is considered. The surface of the poroelastic wspace

is indented by a rigid circular indentor of diameter

20 which applies a prescribed total load of Po with tirne variation in form of a Heaviside step function as follows:

(a) Cylindncal crack

(b) Pemy-shaped crack Figure 4.6 Surface indentation of a cracked poroelastic half-space.

where

-

Figure 4.7 Finite elernent discretkation of poroelastic half-space.

is the normal contact stress between the geomaterial and the rigid indentor. The

region contains a cylindrical crack of length L located at a finite distance b h m the center of the rigid indentor (Figure 4.6a). The faces of the fiacture are considend to be smooth and permeable. Figure 4.7 illustrates the typical finite element discrethion used in the

numerical rnodeling associated with the fhcture d e k e d by the geometricai parameters

Uu=2 and b/a= 1. The location of the outer boundary of the poroelastic medium is specined in relation to the dimensions of the l d e d region. Suitable mesh refinements are incofporated in the

vicinity of the crack tip and in the loaded region in order to improve the accuracy of computational estimates. The soi1 skeleton and pore fluid are assumed to be incompressible. The results of particular geotechnical interest relate to the effect of the singularity at the

crack tip and the geomeincal charactenstics of the crack on the tirne-dependent degree of

which is based on senlement of consolidation settlement. The degree of consolidation (0 rigid indentor is defined as follows

where

uO(t) = u.(O,O,r).

The tirne factor associated with the consolidation process is

selected as

where

Figure 4.8 illustrates the influence of the crack tip stress singularity and pore pressure boundary conditions of crack faces on the degree of consolidation and the pore pressure

filed at the crack tip correspondhg to geometrical parameten Uu=l and b b l . The results indicate that the stress singularity of the crack tip does not greatly influence the resuits, whereas the pore pressure boundary condition of the crack faces has significant

influence on the local and global poroelastic bebaviour.

4

with singularity

E

Y

Ca

0 0.6 -

0.001

0.01

I 10 T (Tirne factor)

0.1

100 1000

Figure 4.8 Effect of stress singuiarhy and crack penneability on (a) consolidation

setîiement of the indentor; and (b) pore pcsninat the crack tip.

Figure 4.9 Effect of geometrical characteristics of the cylhdrical crack on the

consoliâation settlement of the indentor.

Figure 4.9 illustrates the variations in the degree of consolidation settlement of the indentor and the manner in which the poroelastic response is influenced by the length and location of the cylindncal crack. Figure 4.9a indicates that the consolidation rate increases with crack length (Ua)for the case of b/o=l. This shows the influence of cracks on alteration of both the defomability and hydraulic conductivity characteristics of the saturated pomus medium. Figure 4.9b indicates that the greatest change in the poroelastic behaviour of rigid indentor is observed when the crack is located right at the edge of contact point (b/cr=l). For cracks located away from the edge of indentor (bh> or 4) the effect of crack location on the poroelastic behaviour decays with distance. 44.2 Penny-Sbaped Crack

This section examines the behaviour of a rigid anchor bonded to the surface of a poroelastic half-space containhg a Bat penny-shaped CO-axial crack located at the depth of h (Figure 4.6b) from the surface. Figure 4.10 illustrates the variation of degree of consolidation for different depth ratios Nb for both impermeable and permeable bonded interfaces. The numerical results are cornpared with the analytical solutions given by Yue

and Selvadurai (1995) for an intact poroelastic material. Influences of crack geometry charactenstics, sirnilar to those obtained in section 4.4.1, are observed for the poroelastic behaviour of a bonded anchor.

- -andAnalp.cal solution Selvadurai Yue

(1995)

/ /J

Intact; UZa=O

(a) Impermeable bonded interface

-

--

Analytical solution Yue and Selvadurai ( 1995)

@) Penneable bonded interface

Figure 4.10 Effect of geometricai characteristics of the penny-shaped crack on the

cullsolidation settlement of the indentor.

CHAPTER 5

TRANSIENT MOMNG BOUNDARY PROBLEMS

IN POROELASTIC M E D U

5.1 Introduction

The theory of poroelasticity has wide applications in the study of a variety of problems in

geomechanics, biomechanics, materials engineering, en\ironmental geomechanics and energy resources recovery fiom geological formations. A survey of recent developments

in the applications of the theory in these areas is presented by Selvadurai (1996). A great majonty of these problems focus purely on initial boundary value problems where the

boundary conditions are kept Lxed both spatially and temporally. The boundary conditions can relate to either the tractions exerted by the porous skeleton, or the displacements of the porous elastic skeleton or the pore fluid pressure and its spatial

denvatives. There are a significant nwnber of problems which specifïcaily relate to moving boundaries where the boundary conditions themselves are time dependent. In this

Chapter we examine the time-dependent moving boundary problems associated with crack extension in proelastic media

In fust setion of the Chapter, we adopt a computational procedure to examine the g e n e d ~ u e d - m o d emoving boundary problems arising fiom the loading of cracks in poroelastic media. It is shown that the fdte element technique can be conveniently adopted to develop incremental approaches for the solution of such moving boundary problems. in particular the computational methodology accounts for singular stress fields at the crack tip. The computational scheme accommodates the situations where the quasistatic crack extension takes place dong a trajectory govemed by a mixed-mode crack extension concept applicable to the porous skeleton of the proelastic medium. The incremental nature of the iterative scheme allows for the tirnedependent analysis of the problem where the pore pressure and tractions in the geomatenal fabnc and displacements are appropriately adjusted. The searching scheme for the identification of orientation of crack extension is achieved through a double-node splitting technique either at inter-elernent boundaries or through an element. Due to the tirne-dependent nature of the problem, the computational procedure cannot be readily incorporated within a mesh adaptivity scheme. For this reason, the quadratic finite elemcnt discretizations

with a pemissible degree of refuiement are used in the computational modelling. The numericd procedure also prevents overlapping at crack boundaries throughout the time-

dependent anaiysis and accounts for the continuai updating of the boundary conditions.

The numerical scheme is used to obtain solutions to a class of problems of interest to geomechanics and the accuracy of the cornputational scheme has ken verified witb known analytical solutions and limiting cases recovered through analogous problems in classicai elasticity.

The tirst category of problems examined here deals with the situations where symmetry conditions imposed a priori dictate the path of the crack extension. It examines the

moving boundary problem arising fiom the symmetric indentation of a penny-shaped crack smoothly embedded in a proelastic medium, by porous rigid circular discs. The analysis focuses on the study of the stable quasi-static crack growth that results fiom the intemal loading of the penny-shaped crack. The self-simiiar crack extension results in the moving of the crack tip and the alteration of boundary conditions goveming the tractions,

displacements and pore fluid pressures in the incremental crack extension zone. The analysis focuses on the evaluation of timedependent crack extension when the disc inclusions are subjected to either constant imposed displacements or co-t

imposed

total loads with variable loading rates (see e.g. Eq. 2.20). The accuracy of the computational results is also verified by appeal to known analytical solutions applicable to limiting cases associated with elastic undrained or elastic fully drained cases.

The second type of problem relates to the examination of non-planar quasi-static extension of a cylindrical crack emanating from the edge of a ngid circular indentor into a poroelastic half-space. The extension of the crack follows a mixed-mode crack extension

criterion applicable to the porous skeleton. The pore pressure, displacement and traction boundary conditions are modified in the crack extension zone. The numerical results

presented for a variety of probIems indicate the influence of crack-indentor interaction on the stress intensity factors and coosoiidation rates and the effect of undrained compressibility of a geomatenal on the crack extension pattern. 5.2 Crack Extension Modelling

The fuiite element method described in Chapter 4 can be extended to examine the crack extension problems in poroelastic media in order to deveiop the cornputational mode1 it is necessary to establish a crack extension criterion applicable to the porous solid skeleton. The subject of fracture extension in a briale elastic solid bas b e n studied very extensively over the past five decades. Such studies have been motivated by the interest in the examination of crack extension in both rnetailic materials such as steel, and geomaterials such as concrete, rock and ice. Extensive accounts of these developments

can be found in the literature on fiacture mechanics (see, cg. LiebowiR 1968; Atkinson, 1979; Broek, 1982; Sih, 1991). In studies related to the crack extension in brittîe elastic materials, it is necessary to postdate two cnteria The crack extension cnterion establishes the stress conditions necessary for the onset of crack extension. The second relates to the cnterion which establishes the orientation of crack growth.

5.2.1 Criteria for Onset of Crack Extension

The onset of crack extension in brinle elastic solids cm be described by variety of cnteria. Such cntena are invariably developed on the basis of experimental investigations on fracture toughness testing of geomatenals such as concrete mortar and rock. A simple fom of a criterion for the onset of crack extension can be expressed in tems of the fracture toughness of the material in the crack opening mode. Accordingly, in these

deveiopments, the crack extension in a poroelastic medium is assumed to take place when the stress intensity factors applicable to the singular effective stresses at the crack tip

satisfies the condition:

where Kit is the cntical value of the stress intensity factor in the crack opening mode.

The result (5.1) can be generalized to inciude the influence of mode II or flaw shearing effects. The simple mixed-mode fiacture criterion for isotropie brinle elastic solids given by Erdogan and Sih ( 1 963) takes the form

where 8 is the orientation of crack extension which can be obtained fiom the crack orientation criterion.

The midies by Sih (1974) indicate that a generalized theory for the onset of crack extension can also be posed in relation to the local strain energy density at the crack tip. The theory does not require the calculation of energy release rate and thus possesses the

inherent advantage of king able to accommodate crack extension processes in which al1 modes of crack extension contribute to the local stmb energy density function. The strain

energy density function S for a state of plane strain at the crack boundary can be written

as S = alJi

+ Za,&K, +a&

Where

+cos8)(3-4"-cos8)

a,,= '(1 1 6 ~

1

a,,= -sin 8[4 cos 8(1- 2v)I 16ci 1 a,,= -[4(1161i

v)(l COS^)+(^ + C O S ~ ) ( ~-Cl)]O S ~

It cm be shown that the stationary value of Sm. can be used as an intriasic material

parameter, the value of which at the onset of crack extension Sc, is independent of the crack geometry and loading.

5.2.2 Orientation of Crack Growth

In general the two stress intensity factors KIand KI,are present at the crack tip location. Consequently, a generalized crack extension critenon should incorporate the infiuence of

both stress intensity factors. In this study the orientation of crack growth is examiaed by

employing the criteria postulated by Erdogan and Sih (1 963).

The maximum stress criterion assumes that the crack will start to extend in the plane which is normal to the maximum circumfrential stress am (Le. stress state refened to the local polar coordinate reference located at the crack tip) in accordance with the condition

for determinhg the initial angle of crack growth 0. This criterion has been used quite extensively for the study of quasi-static crack extension paths in bririle elastic matenals. Recently, ten Busschen and Selvadurai (1995) and Selvadurai and ten Busschen (1995)

have applied these criteria to the study of hgmentation of fibres embedded in brittle eiastic solids and the extension of cracks from the f'ragmented fibre locations. Excellent agreement is obtained between computational simulation of crack extension paths and

experimental results. Further developments of crack orientation criteria for the crack extension in mixed mode are given by Sih and Theocaris(1 979) and Sih (199 1).

5.2.3 Crack Modelling Techniques

Extensive research effort into modelling the crack extension problem has resulted in the development of various cracking rnodels which cm be generally classified in two categories of discrete crack models and smeared crack models. These models can be applied either with a strength-based or tiacture mechanics based crack extension criterion.

In the discrete crack model, the cracks are simulated by separation of nodal points in finite element mesh. This method is generally implemented by disconnecting the nodal points for the adjoining elements (Figure 5.1) which results in a discontinuous displacement field at crack boundaries. This procedure was fm introduced by Ngo and Scordelis (1967) and then modified by Nilson (1968) in order to allow the finite element model to generate the location of cracks. ln these developments, the discrete crack is modeled by utilizing two nodes at one point (double nodes) dong a predefined crack path

in the finite element discretizatioa. The crack can extend dong the mesh lines by splitting the double nodes. These developments were related to the strength-based crack extension criteria. Shih et al. (1976, 1979) adopted a sirnilar approach fot the fracture-based extension criteria by shifting the crack tip node in the direction of crack growth.

discrete cracks

cracks

Figure 5.1 Idealization of cracks; (a) nodal separaiion using two or four coincident

nodes; (b) smeared crack (ASCE, 1982).

The main difficulty with the discrete crack model is that as the crack extends, the finite element mesh topology has to be modified to reflect the new configuration of the crack particularly when the crack extension path is not known a priori. This can be approached

by either local remeshing near the crack tip or remeshing the entire domain at each crack extension increment. Accounts of development in these methods are given by Ingdea (1977% b), Shephard et al. (1985), Vdliappan and Murti (1985), and Wawryaek and

IngrafTea (1 989) in connection with the FRANC cornputer program.

The smeared crack model fmt introduced by Rashid (1968) assumes that the cracks are smeared in a continuous manner with an infinite number of parallel cracks within the

finite element (Figure 5.1). There is no need to redefine the fmite element mesh topology afier extension of the crack. This method has been widely used in the finite element analysis of concrete structures (ASCE, 1991). This mode1 is based on a concept similar to that employed in the description of the theory of continuum damage mechanics

(Kachanov, 1958) discussed in Chapter 7.

In this study, we adopt and extend the discrete crack model for elastic materials given by Shih

et

al. (1976, 1979) to examine the quasi-static mixed-mode crack extension

problems in poroelastic media This approach allows modelling the pore pressure boundary conditions for the extending cracks more explicitly and realisticaily.

Application of remeshing techniques for this class of transient problems is computing

intensive and relatively inaccurate since dl the idormation at nodal points ffom the previous time step shodd be mapped to the updated finite element discretization. The node splitting technique with the modification of dowiog the elements to be subdivided into two separate elements is employed in the computational scheme.

53 A Computationai Mode1 for Crack Extension ï b e computational methodologies that are developed in Chapter 4 for the study of stable

cracks can be extended to examine the tirnedependent problem of the mixed-modecrack

extension in poroelastic media. The incrementai nature of the iterative scherne ailows for the time-dependent anaiysis of the problem where the pore pressure and tractions in the geomaterial fabnc and displacements are appropriately adjusted. 'Ihe quasi-static crack extension is simulated by introducing the displacernent discontinuity employing the double-node technique (Nikon, 1968). This technique requires that in zones which are prone to cracking, each node in the finite elernent discretization be implemented with double nodes having the same coordinates. During the anaiysis whenever the crack extension criteria are satisfied, the double nodes split into two separate nodes resulting in a displacement discontinuity. In this technique, the

dimension of the general stiffiess matrix will not change as a result of crack extension. This aspect enhances the computational efficiency.

The stress intensity factors related to modes 1 and II are calculated by employing rhe displacernent correlation method (Equation 4.4). The mixed-mode crack extension criteria given by Erdogan and Sih (1963) are assumed to govem the crack growth in the porous fabnc (Equations 5.2 & 5.4). If these conditions are satisfied then, the crack tip moves to the new location in the finite element discretization which remains unchanged

The approximate location of the crack tip is obtained by either splitting double nodes ar inter-element boundaries or subdividing the crack tip element (Figure 5.2) depending upon the crack orientation. The mid-side nodes close to the crack tip are then placed at the quarter point to generate the new crack tip elements. A similar searching scheme for the identification of crack orientation has been used by I n w e a (1 977% b) in comection with the elastic analysis of crack extension. This iterative procedure continues in an incremental fashion until a stable state of a crack is achieved at each time step of the analysis.

The displacements, tractions, and pore pressure boundary conditions change depending upon the extent of crack opening. For the opening mode (KI>O) of crack extension, the pore pressure boundary condition is adjusted in the opened region to reflect a zen,pore

Figure 5.2 Identification of crack orientation.

rigid links with a unilateral constraint

Figure 5.3 Interaction conditions on the faces of the crack.

(i) Cornpute Kiand KI,(Eqs. 4.4a & b): IF K,< O ; Add rigid links between crack faces. Retum to (i).

ELSE;

Proceed to (ii).

(ii) Check fracture criterion:

-

3Kll cos- - - 2 2K,, 78

-

K,sin 0 - KI, (3 cos û - 1) = O

(Eq.5.4)

Yes: Crack extension. Proceed to (iii). No: No M e r crack extension. Exit.

(iii) Locate new crack tip (CT) by following approximations: X

IF O 5 8 < - ; SET 0 = O 8

X

3X

8

8

IF - 1 0 < -;

= CT moves to node k (Fig. 5.2).

SET O = - = CTmoves tonodem. X

4

Split double nodes and generate quarter point singular elements. (iv) Update pore pressure boundary conditions at crack faces:

IF Ki > O IF

K'=O

SET p= O (permeable). SET

PA

= PE

PB= P D

(Fig. 5.2)

continuous pore pressure pressure

Box 5.1 Crack extension algorithm.

pressure boundary condition (permeable crack). Altematively, it is assumed that the pore pressure field remains continuous for cracks which extend in only shear mode (K,=O,

K,,+O) but remain closed (impermeable crack). Al1 the cornputations are perfomed for each time step of the time-dependent analysis.

The interaction conditions are incorporated on the faces of the crack to prevent overlapping of the crack boundaries. This is achieved by using rigid links with a unilateral constra.intcrossing the faces of crack in case of overlapping (K,4)throughout the time-dependent analysis (Figure 5.3). The rigid links disappear when the interaction forces are tensile. The details of the crack extension algorithm are given in Box 5.1. The computational scheme is used to obtain solutions to

MO

problems of interest to

geomechanics and the accuracy of the numerical scheme has k e n verified with known

analytical solutions and Iimiting cases recovered through analogous problerns in classicai elasticity. 5.4 Intemal Indentation of a PennyShaped Crack

The problern of the opening of a penny-shaped crack in a saturated poroelastic geomaterial by smooth porous disc inclusions indenting the faces of the crack is

examined in this section. The class of elasticity problems, which examines the traction loading of a penny-shaped crack in an elastic solid, was examined in detail by Sneddon

(1946) and Barenblatt (1962) (see also Kassir and Sih, 1975; and Cherepanov, 1979). Selvadurai and Singh (1984) and Selvadurai (1985) examined the associated elasticity problem where the surfaces of the crack were subjected to displacement constraints. The motivation for the problem stems fiom its applicability of such a mode1 to the study of

the opening of hydraulic hctures in resources bearing geological media. Selvadurai (1994% b) bas recently extended these studies to include pre-hctured bi-material

geological media (see also GladweU, 1995).

Figure 5.4 Intemal indentation of a penny-shaped crack. We consider the axisymmetrk problern related to the symmetric indentation of a penny-

shaped crack, of radius b, located in an infinite poroelastic medium by a smoothiy embedded porous penny-shaped inclusions of radius a (Figure 5.4). The symmetry o f the

problem about the plane t 0 allows us to consider a single half-space region occupying 1 >O. The rigid inclusion is subjected to either a prescnbed total load or a prescribed

displacement with a time variation in fom of a Heaviside step function and a ramp b c t i o n loading as follows

where H(t) and G(t) are the tirne fiinctions (Figure 5.4).

The application of the symmeaic loading to the inclusion is assumed to initiate self similar extension of the crack when the stress intensity factor applicable to the soi1

skeletal stresses analns a critical value. During this self similar expansion, the

details at A permeable

Figure 5.5 Finite element discretization of indentation problem. displacement, traction and pore fluid pressure boundary conditions change and this results in a movhg boundary problem. The analysis focuses on the study of the quasi-static

crack growth that results fiom the internai indentation of the penny-shaped crack. The

extension of the crack is govemed by the anainment of a critical stress intensity factor at the crack tip (Equation 5.1). The crack extension results in the moving of the crack tip and the alteration of boundary conditions governing the tractions, displacements and pore

fluid pressures in the incremental crack extension zone. The analysis focuses on the calculation of the tirne-dependent evolution of crack extension. Figure 5.5 illustrates the

finite element discretization used in the numerical modelling. The soi1 skeleton and pore fluid of the porous medium are assumed to be incompressible (vu+ 0.5). The result of primary interest to the crack extension modelling relates to the evaluation of

the time variation of the stress intensiîy factor for a stationary crack problem and the time

Figure 5.6 Variations of (a) the crack opening stress intensity factor;

and (b) evolution of crack length for a prescribed total load P.

Figure 5.7 Time-dependent Ioad relaxation and stress intensity factor for

a stationary crack subjected to a prescribed displacement.

Figure 5.8 Tirnedependent load relaxation and thnedependent evolution

of crack length for a prescribed displacernent loading.

variation of crack length for a quasi-static crack extension problem. Figures 5.6 illustrates the time variation of stress intensity factor K,for a Heaviside type of loading. The upper

asyrnptotic value corresponding to the drained elastic behaviour shown in Figure 5.6a is derived fiom the analytical solution given by Selvadurai (1994a). The lower asymptotic

value corresponding to the undrained elastic behaviour is given by Atkinson and Ciaster (1 99 1) as follows

It is s h o w that the resuits for the poroelastic behaviour of discs are consistent

\cith

Iimiting analytical solutions. In particular, it is observed that the transition fiom r ~-WO

the

d to

need not be rnonotonic. The influence of pore fluid migration can result in the

alteration of the stress intensity factors at the crack tip in a poroelastic medium. The fluid migration towards the crack tip leads to a volume expansion of material which tends to initially close the crack in the vicinity of crack tip. This results in a decrease in the stress intensity factor. The crack, however, begins to open up as the pore pressure ciiffision takes place into the permeable faces of the crack. Figure 5.6b presents the variation of the crack length evolution for the situation where the crack is allowed to experience self similar extension. In these circumstances, it is

observed that crack extension in to stable states (i.e. locations where K,/K,,
Figure 5.9 Surface indentation of a poroelastic half-space.

and a permeability of k=2x10*'~ mis, e.g. granite). Figures 5.8 illustrates the computationai results related to an extending crack. The innuence of the rate of loading on the load relaxation and the final crack length cm be observed. For slow rates of

loading, the final length of the crack tends to be smaller. 5.5 Indeata tion of a Poroelastic Half-Space

In this section the computational modelling of the axiqmmetric problem related to the surface indentation of a poroelastic half-space region by a rigid circular indentor is examined. The non-planar quasi-static extension of a cylindncal crack emanating from region at the edge of the rigid indentor into the poroelastic half-space is examined. The extension of the crack follows a rnixed mode crack extension cnterion (Equations 52 &

5.4) applicable to the porous skeleton. The pore pressure. displacement and traction boundary conditions are rnodified in the crack extension zone. The numerical resuiu presented examine the influence of the crack extension on the consolidation rates and influence of undrained compressibility of the geomaterial on the crack extension pattern. The surface of the poroelastic half-space is indented by a miooth rigid indentor of radius a (Figure 5.9)

-

t-igid disc \

r permeable

a

impermeable

1

-

F

*

1Ocr

Figure 5.10 Finite element discretization. which is subjected to a load Po) that is kept constant (i.e. a Heaviside step function of total load). The faces of the developed fracture are considered to be smooth and

permeable. The location of the outer boundary of the proelastic medium is specified in relation to the radius of the rigid indentor. Figure 5.10 illustrates the fuiite element discretization and associated boundary conditions used in the numerical modeling. The soi1 skeleton and the pore fluid are assurned to be compressible with parameters v=0.2, ~ " 4 . 3A. critical stress intensity

(a) Mesh A

(b) Mesh B

Figure 5.1 1 Crack extension patterns for two mesh discretizations A and B.

Figure 5.12 Effict of crack extension on the degree of consolidation.

Figure 5.13 Effect of undrained compressibility on the crack extension path.

factor KIc=l.O MP&

is assurned for the geornaterial fabric. A cylindncal starter crack

with length of O.la at the edge of the indentor is considered. The incremental starter

cracks with either basic conical or penny-shaped configurations are used by Selvadurai and ten Busschen (1995) for the elastic analysis of crack extension with axial symnetry.

Alternative procedures require the assessment of the initiation of a shear crack by appeal to the stress analysis of an intact solid subjected to indentation (Selvadurai, 1997a). The

result of primary interest relates to the assessment of the influence of the crack extension on the time-dependent degree of consolidation senlement. The degree of consolidation

(0which is based on sealement of rigid indentor can be defined as follows

where

up(t) = +(O,O,r).

The finite element dixretization in the vicinity of the edge of

the indentor is refmed to examine the mesh objectivity of the results. Figure 5.1 1

illustrates the crack extension patterns associated with two rnesh contigurations .4 and B,

which are found to be similar. However, the geornetries of crack extension fiom these numencai experiments are slightly different which could be attributed to numerical procedures used to evaiuate crack extension and orientation of crack extension. The timedependent consolidation settiement of the indentor shown in Figure 5.12, indicates that the effect of crack extension on the consolidation rate is nominal. The numerical analysis

shows that the entire crack extension process ?&es place instantaneously and M e r pore pressure difision does not change the crack extension geometry.

The effect of undrained compressibility of the geornaterial on the orientation of crack path is illustrated in Figure 5.13. The crack tends to penetrate more into the porous

medium for more compressible rnaterials. For an hcompressible material

(vu-, OS),

there is almost no crack extension developing in the porous medium. These resuits agree

with analytical/compuiationd solutions given by Selvadurai (1997a) for brittle crack extension fiom the boundary of a circular punch indenting an isotropie elastic half-space.

CHAPTER 6

STEADY MOVING BOUNDARY PROBLEMS

IN PORUELASTIC MEDIA

6.1 Introduction

The computational modelling of hcture problems in poroelastic media can be extended to examine the steady staie problem of a crack propagating through an infinite saturated

porous medium at a finite velocity. Mien the assumptions of steady state crack extension are invoked the governing equations of poroelasticity are modified. The modification

talces the form of the elimination of the time variable by a suitable transformation which

accounts for the propagation velwity at the crack tip. In this Chapter, we present the basic equations goverring the steady state propagation of a semi-infinite crack under conditions of plane strain and axial symmetry.

The problem of steady state constant velocity crack growth in an isotropie elastic medium was first examineci by Yoffe (195 1) for the case of plane strain deformation. YofTe ( 1 95 1)

developed the analytical remlts for a crack of fixed length propagating in an elastic

medium, which is subjected to uniforni remote tensile stresses (Figure 6.1), fiom the solution of the goveming elastodynamic equations. Radok (1956) and Broberg (1960)

extended this study to examine the steady self-similu extension of a crack in an elastic material. The limits of propagation velocity in these studies are established in relation to Rayleigh wave velocity in elastic materials. A recent review of these developments is presented by Freund (1990). The classical fracture mechanics considerations for elastic materials indicate that when the energy release rate of the fracture is greater than that of dnving mechanism, the

extension of the crack is unstable and the failure is rapid. In general the initiation and extension of cracks in elastic materials are transient phenomena. However, a steady state of crack extension can be obtained for lirniting times when the propagation velocity at the

crack tip becomes constant. To preserve the conditions of steady state of crack propagation, the confiiguration of an elastic medium and the üaction distribution should be time invariant in a reference coordinate system which moves with the crack tip. For

poroelastic rnaterials, the flow of energy into the pore fluid tends to stabilize the crack growth and that can result in a quasi-static extension of the crack with a certain velocity.

The phenomena of quasi-static crack propagation in poroelastic geornaterials cm be

atîributed to development of landslides in overconsolidated clay (Palmer and Rice, 1973; and Rice and Cleary, 1976) and aftenhock events in an earthquake (Booker, 1974). The

quasi-static crack growth also govem the hydradic fiacturing phenomena used quite extensively in oil resources recovery (Ruina, 1978; Huang and Russell, 1%!Sa, b; Boone and Ingraffea, 1988). Rice and Simons (1976) and Simons (1977)have given analytical

solutions to problems of a semi-infinite crack propagating quasi-statically in shear mode through a porous medium for various types of pore pressure boundary conditions on the crack faces. Cheng and Liggen (1984b)applied the boundary integral equation method to

solve a sirnilar problern. Recently Craster and Atkinson (1991,1996) have examined in more detail the crack tip stress and pore pressure fields for steadily propagating semiinfinite poroelastic cracks. They have given the analytical solutions for the pore pressure and stress fields near the crack tip for various boundary conditions that are applied to the

pore pressures and the tractions on the crack faces. While the class of two-dimensional

plane strain problems give nse to steady state crack extension phenornena, the equivaient ciass of problems involving extension of circular cracks do not ykld a steady state.

In this Chapter the fuiite element formulation of the transformed quatiom governing the steady state propagation of a crack in a poroelastic medium is first developed by employing the Galerkin technique. The f i t e element approximation results in a system

of non-synimeîric coupled maîrix equotîom which are velocitydependent. The cornputational scheme is then verified by appeal to analytical solutions for the pore

pressure and displacement fields at the crack tip given by Rice and Simons (1976) and Craster and Atkinson (1991). The versatility of the computational modelling procedure is demonstrated by the application of the methodology to the study of the opening of a plane crack by a dipole of rnoving forces or by a smooth rigid punch of uniforni thickness. It is s h o w that the stress intensity factors at the crack tip due to the steady state movement of

the crack tip can be evaluated quite conveniently by the numerical scheme. The

capabilities of the computational procedure are M e r established by considering, for the fiat time, a crack extension associated with a problem exhibiting axial symmetry. This

involves the penetration of a thin cylindrical shell with a *blunt &ont" into a brinle saturated geomaterial. The computational procedure gives estimates for the stress intensity factor Kiwhich can be cornpared with analogous resdts for the rwo-dimensional plane strain problem.

6.2 Steady Crack Propagation

An examination of the Iiterature on elastodynamic fracture mechanics indicate that the first mathematical treatment of steady crack propagation under plane strain conditions was conducted by Yoffe (1951). She examined the problem of a nraight crack of fixed

length 2a propagating through an infinite elastic medium at a constant velocity V (Figure

6.1). The medium is subjected to a remote uniform tension field a,. Yoffe(1951), through solution of the elastodynamic governing equations, obsen-edthat if the crack

Figure 6.1 Steady crack growth problem examined by Yoffe (195 1). propagates in a direction normal to the maximum tensiie stress, there is a criticai at which the crack tends to become c w e d or propagation velocity of about 0.6~~

branched out. We note that c~ is the velocity of transverse (or shear) waves propagating through an elastic body which is given by

where p is the mass density and p is the shear modulus of an elastic material.

The scope of mathematical modelling of dynamic fracture phenornena was then extended by Broberg (1960) to examine the solution to the steady problem of self-similar

expansion of a crack in a uniforni tension stress field in an elastic medium. Based on the

positive definiteness of the energy flow to the crack tip, Broberg (1 964, 1989) proved that cracks cannot propagate in crack opening mode 1 with vefocities higher than the Rayleigh velocity CR given by

where a. is a constant which depends on the value of Poisson's ratio. This constant is

constrained to range in a narrow domain close to unity (e.g. ao=0.88for v=O and ao=0.96 for v=O.5) and c m then be assumed equal to unity for most engineering applications (Davis and Selvadurai, 1996). Most of analytical solutions for the steady crack growth in poroelastic media appear to bave negiected these Iimiting bounds on the propagation

velocity which are attributed to matenal phenomena. It is observed that the ratio of this critical velocity to the permeability coefficient of porous material is a key parameter which characterises the relevant limits for the steady behaviour of crack propagation in

poroelastic media The relevance of these limiting bounds on poroelastic effects of crack propagation shouid be considered. We postulate that a crack located in a poroelastic medium becomes unstable and staris to

propagate if the stress intensity factor KIapplicable io the prous fabric reaches a critical value of KIc. In poroelastic and other materials which exhibit dissipative phenomena the propagation of a ftacture is likely dynamic and unsteady in the initial stages. However, a steady state of crack extension cari be obtained at limiting times when the crack extension has occurred over a long period of t h e . In addition, it is observed that Lhe critical stress

intensity factor (43for the steady state crack extension in proelastic matenals depends on the propagation velocity at the crack tip (see e.g. Rice and Simons, 1976; Simons,

1977; and Craster and Atkinson, 1991). For poroelastic materids, the stress intensity

factor decreases as propagation velocity increases. Therefore it is possible that the crack propagates in a quasi-static fahion with constant velocity such that m e r acceleration would r e d t in lowering the stress intensity factor below the criticai b i t .

The problem of the steady quasi-static propagation of a crack in poroelastic media moving with certain velocity and dnven by tende tractions (Figure 6.2) is examined in

a

Infmite poroelastic medium

crack tip (singular ai,) x

;

Figure 6.2 Steady propagation o f a crack in a proelastic medium.

this section. This steady state crack extension behaviour of saturated materiais is assumed to depend on the velocity of crack propagation. The response of material is drained at low

propagation velocities and undrained at hi@ velocity limits. As a result, the stress field

near the crack tip will be a function of propagation velocity.

6 3 Govemhg Equations The idealized problem of a steady semi-ia6nite plane crack moving in a poroelastic medium is considered. The crack moves along the x-direction with a constant velocity V. We further assume that the coordinate -stem is located at the tip of the moving crack (Figure 6.2). The problern is assumed to be quasi-static and inertial effects and body forces in the medium are neglected. The acceleration of the system is asswned to be zero.

For the description of the problem, we select a crack orientation where the crack is

located dong the x-ais, and the crack is also assumed to move along the x-axis. We consider the transformation X=X-Vt

;

r=Z

(6.3)

where (x, z ) is a coordinate reference system rnoving with the crack tip. The explicit timedependency can be removed by writing

Substituthg Equation 6.4 to governing equations 2.9 results in the following timeindependent equations governing poroelasticity as:

Since time dependency is elirninated tbugh transformation (6.3), for a well posed

boundary value problem, only boundary conditions need to be specified on the variables u,q

and p.

6.4 Finite Elemeat Formulations

The Galerkin approximation technique cm be applied to the goveming equation (6.5b)to transfomi the pariid digerential equation into a discretized rnatrix form. The approximation w d for the displacements u, and pore pressure p are given by Equations 2.17 and 2.18 as follows

where

N", fl correspond to

the nodal shape functions for displacements and pore

pressures. The Galerkin technique is applied to the goveming equation (6.5b) which results in the foilowing weak (weighted residud) form of the equation (6.5b):

where R is the domain of interest, and we note that x,=x; x2=z; and by virtue of the specified orientation of the crack x,=x. Application of Green's theorem to the above equation results in the following:

where B is the boundary of domain R. By substituting the interpolation functions given by Equations 6.6 in the above equation, one obtains:

The above equation c m be written in matruc form as follows:

([W- [ W ) ( u } + ([Hl - [EEI)lp) = IF, 1 where

(6.1 la)

where the compressibility parameten a and P are given in Equations 2.3, and IFq)is the

outward fluid flux through the bomdary B.

The finite element fomulation of the steadily propagating fÎacture in poroelastic media can then be written by cornbining Equation 6.10 with its equivalent equilibrium equation

(Equation 2.19) given by

impermeable (ap/dn=O) or permeable @=0)

--,'#

-

7

permeable @=O)

Figure 6.3 Finite element discretization of steady mck extension. The non-symmetric matrix fonn of discretized equations takes the fom

[ , S C

'

R-EE]{"}=Pl g

where

K = stifiess matrix of the soil sketeton; C = coupling matrix related to interaction between soil and pore fluid; CC = modified coupling matrix;

CB = coupling rnatnx associated with the boundary conditions; EE = modified compressibility matrix of fluid; H = permeabifity matrix;

F = force vectors due to extemal tractions, body forces and flow field; u, p = nodal displacements and pore pressures. 6.5 Verification Exercises

The computational scheme developed for the steady propagating cracks in poroelastic media is calibrated by cornparison with known analytical and numerical solutions. The plane strain problem of a semi-infinite crack embedded in a poroelastic medium is considered. The faces of the crack are first ssubjected to a uniform nomai total stress over a finite distance L from the crack tip. The finite element discretization of the problem and associated boundary conditions are shown in Figure 6.3. The crack tip and associated loading are required to move with a steady velocity Y dong the x-direction. Rice and Simons (1976) and Cheng and Liggett (1984b) have given the analytical and numerical solutions for the variation of the energy release rate G given by following:

It can be noted that the velocity dependence of G will materialize through the calculation

of

4.The value of G falls between two limiting cases of a drained value (Ge)and an

undrained value. The crack extension criterion is satisfied when the energy release rate

-

Analytical, Rice & Simons

Figure 6.4 Crack propagation cntenon for unifonn loading of crack.

reaches a cntical value G, which is dependent on the velocity V. Figure 6.4 illustrates

variation of the normalized cntical energy release rate Ga with the propagation velocity

V. The dimensiodess parameter GJG, approaches unity as V+O and asymptotically

reaches a value i12 (i.e. q=(l-vJ/()/(l-v))as V+-.

The results are compared for the case of

q 4 . 3 3 3 which corresponds to an overconsolidated clay (Cheng and Liggen, 1984b). The

maximum discrepancy of 5% occurs for very hi& velocity where the poroelastic effects become highly localized and pore pressure field shows spatial oscillations. Either a very

fine mesh discretization or a special crack tip element which captures the pore pressure

field applicable to high velocities, in an analytic manner, is needed to address this effect.

In the second problern examined, a semi-infinite crack embedded in an Uifuiite

poroelastic medium is considered. The faces of the crack are subjected to an exponentially decaying normal total stress given by

where H(t) is the Heaviside step function, and a is a constant (Figure 6.3). Atkinson and Craster (1991) have given the closed form solutions for the problern. The variation of flaw opening stress intensity factor Kiwith the velocity V takes the form:

where is the elastic stress intensity factor given by (2a)lRoo; AaJ is a function of nomalized velocity a,=aV/C given by the following:

(i) for a permeable crack

(ii) for an imperneable cmck

Penneable (Analytical) Penneable ('Numerical) Irnpenneable (halytical) Impermeable (Numerical)

Figure 6.5 Variation of the stress intensity factor Ki.

- Penneable (Anaifical) --

Permeable (Numerical) hperrneable (Analytical)

Figure 6.6 Pore pressure distribution ahead of crack tip for velocity u,=1.

Also the generdized consolidation coefficient C is given by

Figure 6.5 illustrates the results for the value of stress intensi~,factor K, derived fiom

analytical solutions and numericd simulations as a function of normalized crack tip velocity aV/C corresponding to material parameters v = 0.3, vu = 0.4. The results indicate that there is g w d agreement between the analytical results and computationai results (maximum discrepancy is 3%). Again in this case, an oscillatory behaviour is observed in the pore pressure field at the crack tip at very bigh velocities. The nominal results also

support Craster and Atkinson's observation that the impermeable pore pressure boundary conditions on crack faces result in a lower stress intensity factor than for the permeable case. This indirectly implies a greater effort is required to extend cracks where the crack faces rernain impermeable. The analytical solutions for the pore pressure fields near the crack tip are aiso given in an explicit form by Atkinson and Craster (199 1). The pore pressure field dong the crack propagation axis (x) ahead of the crack tip takes the following forms:

(i) for a permeable crack

(ii) for an impermeable crack

where

and erfc(x) is the complementary error function given by

Figure 6.6 illustrates the pore pressure distribution ahead of the crack tip for a similar

problem with panmeters v = 0.2 and v, = 0.3, corresponding to a normalized velocity a,=l and for the permeable and impermeable pore pressure boundary conditions on crack

faces. There is good agreement between anaiyticd solutions and numerical results. The results indicate that a significant pore fluid suction field can be developed ahead of the crack tip which reduces the effective stresses in the crack tip region particularly for

permeable cracks. 6.6 Wedghg of a Poroelastic Crack

The class of steady state self-simi1a.r expansion of cracks in poroelastic media is of

particular interest to geotechnical engineering and energy resources recovery fiom geological formations. The numericd procedure is utilized to examine the problem of crack grow* due to symmetric wedging of a crack by ngid indentors in saturated geomaterials under conditions of plane straia and axial syrnmetry.

in silu stress 6,

Figure 6.7 Wedging of a poroelastic medium by rigid indentors.

The rnathematical treatment of steady-state self-similar crack extension under plane strain conditions was first examined by Radok (1956) in solution of a moving punch problem in an infinite elastic medium. Broberg (1964, 1975, 1989) has examined the near-tip fields and directional stability of such crack propagation in elastic material at high velocities

both fiom the point of view of experimental observations and analytical approaches.

e method using quarter-point singular hlelin (1 99 1) has recently used the f ~ t element elements for the evaluation of stability cnteria of wedging in an elastic material. 6.6.1 Phne Strain Moving Punch

The plane strain problem of a steadily moving rigid punch wedging a semi-infite crack

in an infinite poroelastic medium is examined in this section. The indentor of finite width

which moves with velocity V results in a steady self-similar extension of the crack. The crack extension cnteria is governed by attainment of a critical stress intensity factor KIc applicable to the porous fabric which depends on the propagation velocity. The shape of the indentor is assumed to be either a cylinder of diameter D (Le. the contact forces are

Figure 6.8 Effect of crack geometry, in situ stresses, and propagation velocity

on the stress intensity factor for a dipole point force wedging.

semi-infinite punch

/-

(a-)

Figure 6.9 Effect of crack and punch geometry on tbe stress intensity

factor KIfor a wedging rigid strip.

idealized as a dipole of point forces) or a strip punch with a thickness of D and a length of a (Figure 6.7). It is assumed that during steady state crack propagation, the crack tip takes its location at a distance L frorn the edge of indentor. The soi1 skeleton and the pore fluid

of porous medium are assumed to be compressible with material parameten v=0.2 and ~ " 4 . 3 ,respectively. Finite elernent discretization given in Figure 6.3 is used for

simulation purposes. The result of primary interest is the behaviour of crack-indentor interaction with variations in the in situ stresses a. and the indentor geometry. Figure 6.8a illustrates the variation of the crack opening stress intensity factor K[with the geometry of the crack and the in situ stresses a. for a rigid smooth cylinder indentor moving mith velocity of DV/C=O.Ol. This indicates that for higher in situ stresses and tougher matcrials, the resulting crack length is smaller. The effect of propagation velocity

on the crack-indentor interaction behaviour is also illustrated in Figure 6.8b in cornparison with the static result for the problem with equivalent geometry (i.e. V=O).

Figure 6.9 illustrates the crack extension cntenon for a moving ngid strip which moves through a saturated porous medium with velocity of DVIC=O.Ol. It is assumed that the indentor always remains in hl1 contact with the porous medium. Similar crack-indentor interaction behaviour can be observed for various indentor geometries. It is also indicated that length of the wedging indentor (a) in relation to crack length (9has no considerable effects on the steady crack propagation behaviour for a/l >3.

6.6.2 Axisjmmetric Penetration of a Rigid Shell

The axisymmetric problem related to the steady penetration of a ngid smooth thin shell through a sanirated porous medium at a constant velocity (Figure 6.10) is also examined. The penetration results in the extension of a cylindncal crack moving ahead of the advancing shell with a "blunt front". This particular problem is of interest to penetration mechanics of well casings and jacking of pipes in stiff saturated overconsolidated soils

rigid

shell

poroelastic medium

9

rigid

shell

/

C

L

Figure 6.10 Rigid cylindrical shell penetrating through a poroelastic medium.

-

impermeable ( d p l a n 4 )

permeable @=O)

4r

Figure 6.1 1 Finite element discretization of the rigid shell penetration.

0.07

1

,--

Plane strain

Figure 6.12 Variation of the stress intensity factor KIwith radius a for a rigid smooth shell penetrating steadily through a poroelastic medium.

and other fluid saturated geological materials.

The generalized problem examined here deals with a rigid shell of wall thickness of D, and radius of a which penetrates with a steady velocity V through a brittle saturated

porous medium. Figure 6.1 1 illustrates the finite element discretization of the probiem and associated boundary conditions. The crack-shell interaction behaviour corresponding

to a normaiized velocity DV/C=O.OI is illustrated in Figure 6.12 for various shell radii a.

The axisymmetric solution approaches the solution for the associated plane strain problem in large radii (say crlD >1000). This result has certain practical ment in that the relative geometric dimensions of the ngid shell (i.e. radius and thickness), which permits

the consideration of a plane strain solution to an axisymmetric problem, can be identified.

CHAPTER 7

CONTINUUM DAMAGE MECHANICS OF POROELASTIC MEDIA

7.1 Introduction

The assumption of linear elastic behaviour of the porous skeleton is a significant limitation in the application of the classical theory of poroelasticity to brittle geomaterials which could exhibit strain-sofiening and elastic stiffness degradation in the constitutive

behaviour of the soi1 skeleton. This non-linear behaviour can be due to development of microcracks and microvoids in the porous fabric which could aiso alter the permeability characteristics of the porous medium (Figure 7.1). The classical theory of continuum damage mechanics (Kachanov, 1958) can be used to model such damage phenomena in porous saturated materials. In this Chapter we present the development of a fùiite element based computational procedure for the modelling of the effects of damage phenomena on

the behaviour of saturated porous media. Classical continuum damage rnechanics simulates the efTect of microcrack developments on the behaviour of materials pnor to the development of macrocracks (i.e. fractures). In

such modelling, damage is interpreted as a reduction in the stifiess of the material due to the generation of microcracks and other microdefects. In this study attention is restncted primarily to brittle eklrric behaviour of material where the associated darnage processes can be described by an isotropic damage model.

Figure 7.1 Typical stress-strain behaviour and permeability evolution in brittle geomaterials. The poroelastic behaviour of fluid saturated geomaterials can be influenced by the

evolution of darnage in the porous skeleton. This notion of continuum damage is relevant to geomaterials such as soft rocks and overconsolidated clays where progressive softening in an elastic sense can occur due to generation of microvoids or microcracks. The occurrence of microcracks and microvoids in the porous skeleton c m have the irnrnediate effects of altering the elastic stiflhess of the porous skeleton and the permeability characteristics of the porous medium. Figure 7.1 illustrates a typical behaviour of brittle fluid saturated geomaterials expenencing soil skeletal damage. The effect of soil skeletal damage on the tirne-dependent consolidation characteristics of saturated geomaterials can be exarnined by incorporating the concept of continuum

damage mechanics into the classical theory of porwlasticity. This can be achieved by representing the stifniess properties of the soil skeleton and the permeability characteristics of porous medium as bctions of the state of damage in the material. Cheng and Dusseault (1993) have developed a finite element procedure which employs

the concept of damage mechanics to examine the poroelastic behaviour of saturated

geomaterials. They proposed an anisobopic darnage cnterion governing the evolution of

elastic parameters based on experimental observations on rocks. However, this damage mode1 neglects the alteration of pemeabilit. charactenstics of the material which can

result fiom the creation of microcracks and microvoids within the fluid saturated poroelastic material.

The developed computational scheme accounts for the modifications in deformability and permeability charactenstics of materiais. An isotropie damage evolution law is employed in the analysis which is characterized by the dependency of damage parameters on distortional strain invariants. Admittedly, the darnage processes are expected to be highiy anisotropic in nature and could be restricted to localized zones. Different phenomenological damage criteria govemhg the evolution of permeability characteristics

are postulated fiom experimental observations on saturated geomaterials. The numericd procedure is utilized to evaluate the extent to which the poroelastic behaviour of a rigid circdar punch indenting a poroelastic half-space can be intluenced by the damage evolution. 7.2 Principles of Damage Mechanics

The concept of continuum damage mechanics (CDM) is generally attributed to Kachanov (1958) in recognition of his initial studies describing the tertiary creep of solids. The

classical damage mechanics simulates the efTect of intemal degradation of materiais prior to the development of macrocracks (Le. hctures). This theory has been widely used to predict the nonlinear response of

varie^ of materials including rnetals, concrete,

composites, ice, fiozen mils and geologica! materials (see e.g. Krajcinovic and Fonseka, 198 1; Simo and Ju, 1987; Selvadurai and Au, 1991; Cheng and Dusseault, 1993; Selvadurai, 1994c; Hu and Selvadurai, 1995; etc.).

From a physical point of view, damage fan be considered as the development of surface discontinuities in the fonn of microcracks or volume discontinuities in the frlnn of microvoids. The process of damage evolution begins fkom the Wgin state of a material

and ends with hcture of the volume element. The nonlinear behaviour of most brittle materials is attributed to the initiation of new microdefects and gro~zhof existing rnicrodefects. This behaviour is modeled by introducing local darnage variables in the analysis. Damage variables reflet average material degradation at the macro-scale normally associated with the classical continuum description. This facilitates the adaptation of the "damage" concept in the theory of elasticity or in any other theory associated with classicd continuum rnechanics. The coupling of elasticity and damage modrls has been investigated by a number of researchers including Sidoroff (1 98O),

Krajcinovic (1 984), Chow and Wang (1987). and Lemaitre and Chaboche (1990). The basic concepts of the theory of continuum damage mechanics are descnbed in the

following sections.

7.2.1 The Damage Variable

A hindamental postulate in the theory of continuum damage mechanics is that the state of

darnage in a material cannot be easily distinguished macroscopically fiom that of the intact state. Therefore ii becomes necessary to seek an intemal variable to quanti& the phenomena which describes the darnage state of the material. At the scale of microcracks, the damage phenomena results in a discontinuous medium. Kachanov (1958) was the first to introduce a continuous variable related to density of such defects. Figure 7.2 shows a representative volume element of damaged materiai which is of a size large enough to contain many defects, and small enough to be considered as a materiai point within the formulation of continuum mechanics (Davis and Selvadurai, 1996). The overall initial cross sectional area Ad is defined by the normal n. When damage occurs, A,, is reduced to the net area by:

2.The darnage variable D, associated with normal II is defined

Figure 7.2 Representative element of vugin and darnage state of matenal.

where

D8=0

D,=D, O
9

corresponds to the virgin state;

Y

a critical value which corresponds to the Fracture of material;

;

characterizes the damage state.

Hypothesis of lsonopy : In general microdefects are oriented and D, is a function of the

vector o. This leads to a damagr variable of tensorial nature (Lemaitre, 1984). In cases where the damage is direction independent (e.g. damage process results in microvoids with a spherical fom or microcracks which have distribution and orientation without a

preferred direction), the damage can be defmed by appeal to a scalar variable D. 7.2.2 The Net Stress Concept

nie introduction of a darnage variable D leads directly to the concept of a net stress which is the stress defined in relation to the net area. For isotropie damage, the net stress tensor an, is related to the stress tensor a, in the undarnaged state by the following

relationship:

virgin material

I .

damaged material 1

I

equivalent vugin rnaterial

Figure 7.3 Hypothesis of strain equivalence. Hypothesis of Strain Equivalence: Lemaitre (1984) introduced the following hypothzsis of strain equivalence: " A damaged volume of materiai under the applied stress a shows

the same strain response as the virgin state subjected to the net stress 6" " (Figure 7.3). This hypothesis ensures that the constitutive laws applicable to the virgin material are also applicable to the damaged material, with the stress in darnage state being replaced by

the net stress. The constitutive equation for the damaged material which exhibits elastic isotropy and isotropie damage can then be written in the form

where

= p ( D ) is the shear modulus applicable to damaged material which is a function

of the damage variable D. Based on the h>pthesis of strain equivdence. the elastic parameter jï can be obtained by scaling down the corresponding parameter of the virgin material by a factor (1-D). Le.

This implies that Poisson's ratio always remains constant; Le., Y = v . It is understood that this is an added consmint when considering three-dimensional States of stress (Ju, 1990).

7.3 Initiation and Evolution of Damage

The gradual degradation in the constitutive properties of the material is as a resuit of continuing growth of either already existing microdefects or the progressive nucleation of new microdefects. The damage variable is therefore an evolving intemal variable which changes fiom its initial value Do (Do=O in virgin state of material) to its critical value D,

which represents the stage at which macrocracks are formed. For a given state of stress. the extent of damage is an intrinsic property of the material which is charactenzed by a damage evolution law. The damage evolution cnteria can either be postulated by appeal to micromechanics or determined by experiments.

Following observation of expenments conducted on rocks, Cheng and Dusseault (1 993) assumed that damage is a function of the shear strain energy and proposed the following damage evolution equation for rocks

where the equivalent shear strain 5, is defined as

where, ,y

is the octahedral shear strain and Ji is the second invariant of deviatoric

stresses (see e.g. Davis and Selvadurai, 1996) and

and q, y are material constants and c is either 1 or O depending on whether the stress state

corresponds to either loading or udoading. The criticai darnage variable D, is associated with the residual strength of material at the final stage of strain-sofiening behaviour

(Figure 7.1).

7.4 Experimental Results with Indication of Damage The nonlinear constitutive behaviour of most brittle geomaterials is due to the initiation

and coalescence of microcracks in the porous fabnc. The influence of such damage on either elastic stifiess (Le. degradation of elastic moduli) or strength (Le. strain softening)

of geomaterials such rocks and concrete has been observed by Cook (1965), Bieniawski et al. (1967) and Spooner and Dougill (1975). The effect of microcrack developments on

permeability characteristics of such saturated geornateriais has also been observed by Zoback and Byerlee (1975), Samaha and Hover (1992), Shiping el al. (1994), and

Kiyarna et al. (1996).

Figure 7.4 shows the experimental results obtained by Spooner and Dougill (1975) for a series of cyclic uniaxial compression tests conducted on concrete. It indicates that there is a reduction in the elastic modulus ( E ) with increase of deformation and a reduction in the strength of the concrete in form of strain sofiening behaviour in post-failure region with plastic deformation. Simo and Ju (1987) have successhilly applied the concept of darnage

mechanics to simulate the experimental results on concrete obtained by Scavuzzo et al. (1983). They developed an isotropie elasto-plastic darnage mode1 by employing a

Figure 7.4 The -miaxial compression behavior of concrete (Spooner and Dougil!, 1975).

Figure 7.5 Variations of volumetric strain and pemeability coefficient as a function of differential stress (Zoback and Byerlee, 1975).

tensorial fom of darnage variables. Also Cheng and Dusseault (1993) developed an anisotropic damage model using a vectorid representation of damage variable based on experimental observations to model the constitutive behaviour of rocks and concrete. m e effect of microcrack development on the pemeability charactenstics of a fluid saturated geomaterial was fust investigated by Zoback and Byerlee (1975) who conducted triaxial tests on granite. Figure 7.5 illustrates the results for variation of the

hydraulic conductivity as a function of deviator stress (q-q) in the triaxial test. It indicates that the permeability coefficient is first reduced slightly due to closure of preexisting microcracks and void channels as a result of elastic deformations and then begins to increase considerably by the growth of existing microcracks and the nucleation of new

microcracks. Shiping et al. (1994) exarnined the permeability evolution of sandstone in series of complete triaxial stress-strain paths. They observed that permeability charactenstics of matenal can increase nearly by one order of magnitude, prior to the peak of stress-strain cuve and can increase up to two orders of magnitude in the strain-

soflening regime where microcracks tend to localize in shear faults. Kiyama et al. (1 996) observed similar results for permeability evolution of granites in triaxial experiments. In

addition to the experimental observations, there has been only limited work in the context of mathematical and computational treatment of pemeability alterations of saturated geomaterials. 7.5 Effect of Damage on Poroelastic Parameters

The effect of soi1 skeletal darnage on the time-dependent poroelastic behaviour of saturated geomaterials can be examined by incorporating the concept of continuum

damage mechanics into the classical theory of poroelasticity. This is achieved by representation of elastic parameters and permeability characteristics of the porous medium as îùnctions of the state of damage in material. The damage evolution criteria

governing the alteration of stifiess and permeability parameters of the porous material

are described in relation to the damage variable in computational procedures. This approach facilitates the adaptation of the continuum darnage concept in the classical theory of poroelasticity applicable to porous continua. 7.5.1 Elastic Parameters

The constitutive parameters applicable to an isotropic elastic material which is experiencing micromechanicai damage can be reprerented as a function of intact elastic properties by invoking the hypothesis of strain equivalence. This relation can be vjnacn in tensorial form by the generalization of Equation 7.3 for a damaged material which exhibits elastic isotropy and isotropic damage as

where CqUis the elasticity matrix applicable to virgin elastic materiais (see e.g. Davis and Selvadurai, 1996). The Poisson's ratio is assumed to be constant. The damage evolution law cari specify the variation of the damage variable

(D) with the state of strain in a

material. The damage evolution law proposed by Cheng and Dusseault (1993) is

employed in this study to mode1 the elastic stiffness degradation eflects in a material. By integration of Equation 7.5, one can obtain the evolution of damage variable as

where Do is the initial value of damage variable correspondhg to the intact state of a material. Figure 7.6 illustrates the variation of the darnage variable D as a function of material parameters q, y and D,.

Figure 7.6 Variations of darnage variable with material parameters.

7.5.2 Perrneability Characteristics

Development of a damage criteria which can account for alterations in the hydraulic conductivity during evolution of damage in saturated geomaterials is necessary for computational modelling of such phenornena in poroelastic media. The literature in the study of coupling between microcrack developments and permeability evolution in

saturated geomaterials is restricted to experimental observations. There has been only limited experimental work in the context of constitutive modelling of permeability characteristics in damaged porous media (Shiping et d , 1994). In this study two phenomenological constitutive models based on the experirnental results observed by Zoback and Byerlee (1975) and Shiping

et

al. (1994) in rocks, are postulated for the

permeability evolution criteria in poroelastic media.

The slight initial reduction in the hydraulic conductivity of porous materials due to dosure of pore voids and pre-existing microcracks in the elastic range prior to the onset

of microcrack developments is neglected in this study. The hydraulic conductivity (k) is

postulated to have either Iinear or quaciratic variations with respect to equivalent shear strain cd given in Equation 7.6 as follows:

where kd is the hydraulic conductivity applicable to damaged matenals and c,, c?, c, and c, are material constants.

7.6 Finite Element Procedures

The concept of continuum damage mechanics is incorporated in the finite element code

developed in connection with this research for computational modelling of poroelastic media. A computationai procedure is developed ushg the damage evolution critcria governing the elastic deformability and permeability characteristics of poroelastic materials. The scalar damage variables are first obtained at nine Gauss points within finite elements fiom Equation 7.8. The constitutive matrix

eqkl and the hydraulic conductivity

kd are then updated at these locations to incorporate the damage evolution. The

discretized governing equations are then solved to obtain the state of strains at each integration point using the updated matrix Ciikl and parameter kd in incrementai analysis. The coupling between the state of strains and state of damage in each time step is solved by an iterative process. The standard Newton-Raphson (see e.g. Smith and Griffiths,

1988) technique is used for the iteration algorithm in analysis. The adopted convergence

criterion is based on the nom of the evolution of damage variable in relation to a specified tolerance o (see Box 7.1). The computational procedures are performed for each time step in time-dependent

anaiysis. The details of the damage algorithm are given in Box 7.1. The initial -te damage (i.e. Do in Equation 7.8) can also be prescribed for a given element.

of

(i) Compute D at Gauss integration points

(Eq. 7-81 (Eqs. 7.6)

(ii) Update the poroelastic parameters

(Eq. 7.7)

(Eq.7.9a) (Eq. 7.9b)

(iii) Solve the goveming equations for u, p, and calculate the strain tensor

(iv) Check convergence criteria

No: Damage growth. R e m to (i).

Yes: No further damage evolution. Exit.

Box 7.1 Damage evolution algorithm.

131

7.7 Indentation Problem of a Damageable Geomaterial

The axisymrnetric indentation of a damage susceptible poroelastic haif-space by a rigid smooth circular indentor is examined in this section. The poroelastic medium is subjected to a total load P with a time variation in form of a Heaviside step function. The soil skeleton and the pore fluid of the geomaterial are assumed to incompressible (v,=O.j). The material parameters employed in the numerical simulations for an isotropie damage variable adopted fiom Cheng and Dusseault (1993) (for uniaxial compression tests on

sandstone) are as follows:

The hydraulic conductivity is assurned to increase one order of magnitude from the initial virgin state to the peak stress (oma=30 MPa). This is consistent with expenmental observations of Shiping et al. (1994) and Kiyama et al. (1 996). Both linear and quadratic variations of the hydraulic conductivity as a function of distortional strain energy, given by Equations 7.9, are considered. The total load is assumed to be P=635 MN. Figure 7.7

illustrates the stress-strain behaviour and the evolution of the damage variable and hydraulic conductivity for the material considered in numerical simulations.

Figure 7.8 shows the finite element discretization of the problem and associated boundary conditions. Figure 7.9 illustrates the numerical results related to the degree of consolidation settlement of the indentor. For the instance where the nonlinear damage mode1 has no alteration in the permeability parameter, the rate of consolidation decreases due to sofiening of the matenal as a result of microcrack developments. The generation of

damage in the soil skeleton results in a higher displacement and pore pressure. The

,a

A

O

2

= 30 Mpa

4

6

Axial strain

E.

8

1O

x l O-'

Axial strain E, x 1O"

30 (c)

Figure 7.7 The stress-strain behaviour; evolution of damage and permeability

characteristics in uniaxial compression of sandstone (adopted h m Cheng and DusseauIf 1993).

indentor

1Ou

permeable @O)

b

permeable

Figure 7.8 Finite element discretization of indentation problem.

darnaged material therefore requires a longer tirne to achieve the sarne degree of consolidation. When the modification of hydraulic conductivity is taken into consideration, the rate of consolidation increases and the excess pore pressures are dissipated at a faster rate (Figure 7.9). The effect of an increase in hydraulic conductivity ovemdes the effect of reduction in elastic stifhess properties on the consolidation behaviour of the material. This difference is more significant when a linear form of evolution law for the hydraulic conductivity is considered, where, the increase in

hydradic conductivity is higher. Figure 7.10 shows the evolution of the extent of the

a

--

Damage (k= linear)

.-..--Damage (k=const.)

0.001

0.01

0.1

1

10

100

1000

T Figure 7.9 Degree of consolidation settlement of indentor.

0.00001 0.0001 0.001

0.01

0.1

1

10

100

T Figure 7.10 Extent of damaged zone (with D > 0.05) in t h e .

1000

Elastic Damage (k=quad.) Damage (k=linear) Darnage (konst.)

r/a=1 da= 1.5 ;i a \

:I t \

;\

: \

\it

i

----..

.* * . - * -

Elastic Damage (k=quad.) Darnage (k=linear) Damage (k=const.)

Figure 7.11 Pore pressure evolution at different depths.

F i g w 7.12 Evolution of damage variable at the edge of indentor.

damaged zone (i.e. where D > 0.05) uith tirne. It indicates that most of damage phenornena takes place alrnost instantaneously and M e r pore pressure difision does

not change considerably the extent of damage.

Figure 7.1 1 illustrates the evolution of pore pressure at two locations, the center of indentor (r/u=O) and the edge of indentor ( r l d .O) corresponding to a depth of r k l . 5 within the poroelastic half-space. nie damage models predict higher excess pore pressures in porous medium which is consistent with observations by Cheng and

Dusseault (1993).

The evolution of darnage variable with time is shown in Figure 7.12 at the edge of indentor (r/u=l) at a depth of z/a=O. 1.

The numerical results illustrate the importance of the effect of penneability characteristics on the consolidation behaviour of porous media when the ifluence of damage effects is

taken into consideration. The influence of damage generation on the fluid transport

characteristics of saturated geomaterials is more attributed to the alteration of the pemeability characteristics than those of elastic properties of the soi1 fabnc.

The finite element discretization was iefmed M e r in the vicinity of the indentor to examine the mesh objrctivity of the results. The numerical resuits did not change

considerably by the mesh refmement.

CHAPTER 8

CONCLUSIONS AND RECOMMEIYDATIONS

In the preceding Chapters, we have presented the finite element modelling of initiation and extension of fracture, and initiation and evoluiion of darnage phenomena in fluid

saturated porous media for both two-dimensional plane strain and axisymmetric problems. In this final Chapter, we shall summarize the main achievements of this

research and present recommendations for funw work. 8.1 Summary and Conclusions

The satuated geological materials with brittle behaviour c m experience non-elastic phenomena in the soi1 skeleton either in the form of continuum micromechanicd damage by development of rnicrodefects, or in the fom of discrete fractures. Development of

such defects in the porous skeleton of sahuated geomaterials will influence the fluid transport characteristics and the poroelastic behaviour of such materials. One of objectives of this research is to develop a methodology for computational modelling of

fracture and damage phenomena in porous media saturated with compressible pore fluids for both two-dimensional plane strain and axisymmetric situations.

A great majority of applications of the theory of poroelasticity focus purely on initial

boundary value problems where the boundary conditions are kept fixed both spatidly and

ternpodly. The boundary conditions can relate to either the tractions exerted to the porous skeleton, or the displacements of the porous elastic skeleton or the pore fluid

transport. There are a significant number of problems which specifically relate to moving boundaries associated with fracture phenornena where the boundary conditions themselves are time dependent. h this thesis, the class of time-dependent moving boundary pmblems associated with either the transient quasi-static crack extension or the steady state crack propagation in poroelastic media is examined. The literature review on the computational modelling of poroelastic media reveals that the numerical treatment of non-elastic bnttle behaviour of the porous fabric (either in

form of continuum darnage or in form of discrete cracks) and the moving boundary problems arising from the transient and steady state extension of cracks in saturated geomaterials have received limited attention. In order to achieve the objectives of the study, we have completed the following phases of the research program: 1. Development of a finite element code for the two-dimensional consolidation analysis

of poroelastic media s a m t e d with a compressible pore fiuid. Eight noded isoparametric elernents are used to represent the intact regions of geological medium. The polynomial shape functions used to describe the variation of pore pressure field are one order lower

than those used for the displacement fields. It is demonstrated that the spatial oscillations in pore pressure results can be minimized with these element formulations. The computational scheme is verified by appeal to various known analytical solutions given

in the literature for consolidation of poroelastic materials. 2. The computational modelling of transient stationary fractures in poroelastic media is examined. The time-dependent local effects at the crack tip for this class of fiacture problems in sanirated materiais is govemed by the interaction between the flow of pore

fluid and the elastic deformation of the soi1 skeleton. The main numericd developments

and related resuits of interest an summarized in the following:

The concept of linear fracture mechanics can be extended to poroelastic media by incorporation of a stress singularity of order fin for the effective stress fields at the crack tip. This is achieved by ernploying the singular traction quarter point crack tip element developed by Henshell and Shaw (1975) and Barsoum (1976) in conventional finite element analysis of elastic problems. The mid-side nodes close to the crack tip in regular isoparametric elernents are shified to their quarter points to generate the singular crack tip elements. In poroelastic materials, the pore pressure field at the crack tip does not eshibit

spatially singular behaviour. However, in undrained behaviour it shows a temporal singularity as r+0'. The pore pressure field at the crack tip is modeled by regular isoparaxnetric elements. The stress intensity factors at the crack tip which characterize the displacement and effective stress fields near the crack tip region are computed by the displacement correlation method. The numerical results for the tirne-dependent evolution of the stress intensity factors are verified by cornparison with knowvn analytical solutions for various pore pressure boundary conditions at the crack faces.

Numer-ical results indicate that the consideration of a stress singularity at the crack tip

has no considerable effect on the global poroelastic responses such as the degree of consolidation. While it has a considerable effect on the local poroelastic effects such

as the pore pressure field at the crack tip. The pore pressure boundary conditions at the crack faces have significant effects on both local and global poroelastic responses. The initial value of the stress intensity factors (as t+0+ ) corresponds to the elastic undrained behaviour with Poisson's ratio vu and their final values (as t - m ) correspond to the elastic fblly drained behaviour with Poisson's ratio v. It is observed

that the transition of the results h m t+0+ to r-

need not be rnonotonic. The fluid

migration towards the crack tip cm lead to a volume expansion of material in the vicinity of the crack tip. This results in a partial closure of a poroelastic crack with a

reductioa in the stress intensity factor at the early tirnes. However, the crack begins to open up as the pore pressure diffusion takes place into the crack tip.

The impermeable pore pressure boundary conditions on the fracture surfaces result in

a lower stress intensity factor than those for the permeable case. This indirectly implies a greater effort is required to extend cracks in impermeable hctures.

3. A general computational algorithm is developed to examine the rnixed-mode quasistatic extension of cracks in poroelastic media. This class of transient moving boundary

problems arises from the loading of the cracks in saturated matends. It is shown that the finite element technique can be conveniently adopted to conduct Uicremental analyses of

such moving boundary problerns. The main features of the computational scheme along

with related results of interest to poroelastic media c m be summarized in following: The computational scheme accommodates the situations where the quasi-static crack extension can take place along a trajectory dictated by a mixed-mode crack extension concept applicable to the prous skeleton of fluid saturated materials. The incremental nature of the iterative scheme allows the time-dependent d y s i s of the problem where the pore pressure and tractions in the geomaterial fabric and displacements are appropriately adjusted. When the crack extension critena are satisfied, the crack tip moves to its new location

in the tirne-invariant finite elemeot discretization. The searching scheme for the identification of orientation of crack extension is achieved either by splitting the double nodes at interslement boundaries or by subdividing the quadrilateral singular element at the crack tip into two triangular singular elements. depending upon the crack orientation. The mid-side nodes close to the crack tip are shifted to their quarter

points to generate the new crack tip elements. The quacirilateral isoparametric finite elements are used with no mesh adaptivity due to the time dependent nature of the problems.

The displacements, tractions, and pore pressure boundary conditions are altered

depeading upon the extent of crack opening. For the opening mode of crack extension, the pore pressure boundary condition is adjusted in the opened region to

reflect a zero pore pressure boundary condition. Altematively, it is assumed that the pore pressure field remains continuous for cracks which extend in oniy shear mode but rem& closed reflecting an impermeable crack.

The interaction conditions are incorporated on the faces of the crack to prevent overlapping of the crack boundaries. This is achieved by using ngid links with a unilateral constraint at nodal points on the overlapping crack boundaries throughout the the-dependent analysis. The numericai scheme is employed to obtain solutions to certain fracture mechanics

problems associated with poroelastic media. The accuracy of the computational scheme bas been verified with known analytical solutions and limiting cases recovered through analogous problem in classical elasticity. It is observed that most of the crack extension in poroelastic media takes place instantaneously upon the application of extemal loads. Funher, the pore pressure drainage pmcess does not considerably change the geometry of the crack path.

The compressibility of the pore fluid and the soi1 skeleton of poroelastic materials is found to govem the extension and orientation of crack path in such materials. niese

observations agree with the results obtained in the classical elasticity. The rate of loading is found to have a considerable effect on the geometry of the crack

extension for poroelastic materials. For fast rate of loading, the crack tends to extend deeper into the poroelastic medium. 4. The computational modelling of poroelastic fracture problems is extended to examine the steady state crack extension in infiaite poroelastic media at a finite velocity. The goveming equations of poroelasticity are modified to replace the time variable with the propagation velocity at the crack tip. The finite element formulation of the steady crack propagation in poroelastic matenals under conditions of plane strain and axial symmetry is derived. The main computational characteristics of the steady crack propagation mode1

dong with related results of interest to poroelastic media can be summarized in the following:

0

n i e finite element formulation of the governing equations of the steady state crack extension in poroelastic media is developed by employing the Galerkin technique. n i e crack extension is assumed to be quasi-static and inertial ef5ects and body forces

in the medium are neglected. The finite element approximation results in a system of

non-symmetric coupled matrix equations which are dependent on the propagation velocity at the crack tip. 0

The computational scheme is verified by appeal to known analytical solutions to the

pore pressure and displacement fields at the crack tip in poroelastic media. Both permeable and impermeable pore pressure boundary conditions at the crack faces are considered in these venfication exercises. The proelastic behaviour of sahirated materials in the steady state crack extension depends on the propagation velocity at the crack tip. The response of materiai is

drained at low propagation velocities and undrained at high velocity limits. The limits of propagation velocity in poroelastic media are established in relation to the elastic Rayleigh wave velocity and the hydraulic conductivity of such media.

The poroelastic effects at the crack tip become highly localized and pore pressure field shows spatial oscillations at very high velocities. A very fme finite element discretization at the crack tip is required io capture the pore pressure field applicable to high velocities.

Similar to the equivaient transient crack extension phenornena, the pore pressure boundary conditions on the crack faces influence the crack tip behaviour for the steady crack extension in poroelastic media. A greater energy is required to extend the

cracks in impermeable fractures. It is observed that a significant pore fluid suction field cm be developed ahead of the crack tip in poroelastic media which reduces the effective stresses in the crack tip ngion particularly for permeable cracks.

The computational scheme c m be used to examine the self-sirrjlar expansion of a crack in an infihite poroelastic medium by wedging rigid indentors. The results

obtained indicate the influence of in situ stresses at the crack tip and the geometry of indentor on the crack-indentor interaction behaviour.

5. The computational modelling of poroelastic media is extended M e r to examine the behaviour of saturated brinle geomaterials which experience micromechanical damage in the porous fabric. The occurrence of microcracks and microvoids in the porous fabric can have the immediate effects of altering the elastic stifiess and the permeability

characteristics of the porous medium. The concept of Continuum Damage Mcchanics can

be incorporated into the classical theory of poroelasticity to evaluate the extent to which the poroelastic effects can be influenced by the darnage evolution in the porous skeleton.

The main features of the finite element procedure and related results of interest to

poroelastic media c m be surnmarized in the following:

r

The elastic properties of the soi1 skeleton and the permeability characteristics of the porous medium are represented as functions of the state of darnage in material. An isotropic damage critenon goveniing the evolution of elastic modulus of the porous fabric is employed in the computational model. It is characterîzed by the dependency

of damage parameters on distortional strain invariant. Various phenomenologicai damage criteria goveming the evolution of pemeability characteristics are postulated based on available experirnental observations on saturated geomaterials. r

The constitutive elastic and permeability parameters are updated at each integration point in the finite element analysis. The couplhg between the state of deformation and the state of damage at each t h e step is solved by an iterative approach in the time-dependent anaiysis. The standard Newton-Raphson technique is used for the iteration aigorithm. The adopted convergence criterion is based on the nom of the

evolution of damage variable in relation to a specified tolerance. The initial state of darnage can aiso be prescribed.

a

The influence of darnage development on the fluid transport characteristics of

saturated geomaterials is more attributed to the alteration of the permeability characteristics rather than to the evolution of elastic properties of the soi1 skeleton. a

It is obsemed that most of the damage process in poroelastic media takes place instantaneoudy upon the application of extemal loads. Further pore pressure changes do not considerably change the extent of damaged zones.

8.2 Recommendations for Future Work

In the preceding sections , we have summdzed the main achievements of this research. The resemh has provided the computational modelling of the initiation and extension of cracks and initiation and extension of damage phenomena in poroelastic media saturated with compressible pore fluids. The finite element procedures are verified by cornparison

with known analytical solutions to various boundary value problerns related to poroelastic media. The computational mode1 has been successfully applied to various fracture and damage mechanics problems associated with poroelastic media. In the ensuing, we

recommend some specific extensions. The computational models developed for fracture and damage phenomena can be

coupled to examine the nodinear behaviour of poroelastic materials at the crack tip region (process zone). The computational methodology c m be extended to examine three-dimensional hcture mechanics analysis of the poroelastic media. Either a mesh adaptive

technique or non-stnictured triangular mesh discretizations can be employed in finite element procedures to examine the quasi-static crack extension in poroelastic media. The anisotropic darnage models can be employed to examine the inherent anisotropy of darnage phenornena in porous media. The darnage critena goveniuig the evolution

of permeability characteristics of a porous medium need to be determined based on

either experimental procedm or micn>mechanicsconsiderations.

A special crack tip element which can capture the pore pressure field applicable to

high velocities in an analytic manner can be incorporated at the crack tip for the steady crack extension problems in poroelastic media.

References Aboustit B.L.,S.H.Advani and J.K. Lee, 1985, 'Variational principles and finite element

simulations for thenno-elastic consolidation', International Journal of Numerical and Anak'ytica~Methods in Geomechanics, 9,4949.

Agbemge L.K.and H. Deresiewicz 1974, 'Chthe indentation of a consolidating halfspace', Isroel Journal of Technology, 12,322-328. Agbemge L.K. and H. Deresiewicz, 1975, 'The consolidation settlement of a circular footing', Israel Journal of Technology, 13,264-269. ASCE Task Committee, 1982, 'Finite element analysis of reinforced concrete', State of

the Art Report, ASCE, New York, USA.

ASCUACI Committee 447, 1991, 'FUiite element analysis of reinforced concrete smictures II', J. Isenberg, American Society of Civil Engineers, New York, USA.

Atkinson C., 1979, 'Stress singularities and hcture mechanics', Applied Mechanics Reviews, ASME, 32, 123-135.

Atkinson C. and R.V. Craster, 1991, 'Plane strain hcnire in poroelastic media'. Proceedings of Royal Sociev, Series A, 434,605633.

Atiuri S.N. (Ed.), 1991, Compu~atiottulMerho& in Mechaflics of Fracture, Vol. 2,

Comput~tional Methoh in Mechanics, J.D. Achenbach (Ed.), North-Holland, The Netherlands.

Banerjee P.K.and R. Butterfield, 1981, 'Transient flow through p r o u elastic media',

Developments in B o u n d a ~Element Methods-2, P.K. Banejee and R.P. Shaw (Eds.), Applied Science Publishers, England, Chapter 2.

Barenblatt G.I., 1962, 'The mathematical theory of equilibrium cracks in brinle fracture', in Ahances in Applied Mechunics, VII, H.L.Dryden and Th. von Karman (Eds.),

Acadernic Press, New York. 55-129. Barsoum R.S., 1976, 'On the use of isoparametric finite elements in linear fracture

mechanics', International Journal/or Numericd Methods in Engineering, IO1 25-3 7.

Bazant 2.P., 1991, Stability of sttuctures: elastic, inelastic. fiocture. and damage theories, Oxford University Press, New York.

Bieniawski Z.T.,H.G. Denkhaus and U.W. Vogler, 1967, 'Failure of hcnired rock', Internatio~lJournal of Rock Mechanics and Mining Science and Geomechanics Abstracts, 6,323-346.

Biot M.A., 1941, 'General theory of three-dimensional consolidation', Journal of Appiied Physics, 12. 155-164.

Biot M.A., 1955, 'Theory of elaslicity and consolidation for a porous anisotropic solid', Journal of Applied Physics, 26, 182- 185. Biot M.A., 1956, 'General solution of the equation of elasticity and consolidation for a porous material',Journol of Applied Mechanics, ASME, 23,9 1-95.

Blandford GE.,A.R. I n e e a and I.A. Liggett, 198 1, 'Two-dimensional stress intensity factor computations using the boundary element method', In~emutionulJournal for Numerical Methoh N, Engineering, 17,387-404.

Booker J.R., 1974, 'Time dependent s t n h following faulting of a porous medium', Journal cf Geophysical Research, 79,2037-2044.

Booker J.R. and C. Sawidou, 1985, 'Consolidation around a point heat source', Internntional Journal of Numerical and Anafytical Methodr in Geomechanics, 9, 173184.

Booker J.R. and J.C.Srnail, 1975, 'An investigation of the stability of numencal solution of Biot's equations of consolidation', International Journal o/&lids and Structures, I l , 907-9 1 7.

Boone T.J. and A.R. Ingraffea, 1990, 'A numerical procedure for simulation of hydraulicaily-driven fracture propagation in proelastic media', International Journal of Numerical and Analytical Methodr in Geomechanics, 14,2 7-47. Boussinesq J., 1885, 'Application des potentials a I'etude de I'equilibre et du mouvement des solides elastique', Gauthier-Villars, Paris, France. Brebbia C.A., J.C.F.Telles and L.C. Wrobel, 1984, B o u n d v Element Techniques, Springer-Verlag, Berlin. Brebbia C.A. and M.H. Aliabadi, 1991, Numericol Fracture Mechnics. Cornputational Mechanics Publications, U.K. Broberg K.B.,1960,' The propagation of a brittie crack', Archivfur Fysik, 18, 159-192. Broberg K.B.,1964,' On the speed of a bnttle crack', Journal of Applied Mechanics, 31, 546-547.

Broberg K.B., 1975. Reports of 7he Research Inrritute of Strength and Fracture of Materiais, Tohoku University, 11, No. 1. 1-27. Broberg K.B.,1989,' The near-tip field at high crack velocities', Internatiortul Journal of Fracture, 39, 1- 13. Broek D., 1982, Elemeniary Engineering Fracture Mechanics, Martinus Nijhoff Publications, The Hague. Cheng H. and M.B.Dusseaulf 1993, 'Deformation and diffusion behaviour in a solid experiencing damage: A continuou damage mode1 and its numerical implementation',

Internarional Journal of Rock Mechunics and Mining Sciences & Geomechanics Abstracts, 30, No. 7, 1323-133 1.

Cheng A.H.D. and J.A. Liggett, 1984a 'Boundary integral equation method for linear

porous elasticity with applications to soi1 consolidation', International Journal of Numerical Method in Engineering 20: 255-278. Cheng A.H.D. and J.A. Liggett, 1984b. 'Boundary integrai equation rnethod for linear porous elasticity with applications to k t u r e propagation', Infernarional Journal uf Numericd Method in Engineering, 20.279-296. Cherepanov G.P., 1979, Mechanics of Brittfe Fracture, R. de Wit and W.C.Cooley (Tram. & Eds.), McGraw-Hill, New Y o k

Chiarella C. and J.R. Booker, 1975, 'The time-senlement behaviour of a ngid die resting on a deep clay layer', Quarterly J o u d of .Mechics and Applied Mathematics, 28,317-

328.

Chow C.H.and J. Wang, 1987, 'An anisotropic theory of elasticity for continuum darnage mechanics', International JowmI of Fracture, 33,3- 16.

Christian J.T.,1968, 'Undrained stress distribution by numerical methods', Journal of

Soi1 Mechnics and Foundations, ASCE, 94, 1 333-1345. Cleary M.P., 1977, 'Fundamental solutions for a fluid-saturated porous solid', InternationalJournal of Solidr and Structures, 13,785-806.

Cook N.G.W.,1965, 'Failure of rock', Iniernntional Journal of Rock Mechanics and Mining Science ami Geomechanics Abstracts, 2,28943. Craster R.V. and C. Atkinson, 1996, 'Theoretical aspects of h c h i r e in porous elastic

media', Mechanics of Poroelastic Media, A.P.S. Selvadurai (Ed.), Kluwer Academic

Publ., A. A. Dordrecht, The Nethertands. C m e T. and R.B. Wilson, 1977, 'Boundary integral equation method for elastic fracture

mechanics', AFSOR-TR-78-0355, 10-1 1. Cryer C.W., 1963,' A cornparison of the three dimensional consolidation theories of Biot

and Terzagh', QuarterIy Journal of Mechanics and Applied Mriihematics, 16,40 1 -4 1 2.

Darcy H.,1856,'Les fontaines publiques de Dijon', Victor Dalmont, Paris. Davis R.O. and A.P.S. Selvadurai, 1996, EloJlic.).

and Geomechanics, Cambridge

University Press, Cambridge, England Desai C .S. and H.J. Siriwardane, 1984, Constitutive Luws for Engineering M i i o l s with Emphasis on Geologic Materials, Prentice-Hall Inc., Englewood Cliffs, N.J.

Detournay E. and A.H.D. Cheng, 1991, 'Plane strain anaiysis of a stationary hydradic fiacture in a proelastic medium', International Jourmf of Soli& and Shuctures, 27? 1645- 1662.

Dominguez J., 1992, 'Boundary element approach for dynamics proelastic problems', Inrernutional Journalfor Numerical Methoris in Engineering, 35,307-324. Erdogan F. and G.C.Sih, 1963, 'On the crack extension in plates under plane loading and transverse shear', Journal of Basic Engineering, ASME, 85,5 19-527. Freund L.B., 1990, Dynomic Fracture Mechanics, Cambridge University Press, U.S.A. Fung Y.C., 1965, Foundarions ofSolid Mechmics, Prentice Hall, Englewood Cliffs, N.J. Gaiin L.A., 1961. Contact Problem in Theory of Elasticify, LN. Sneddon and H. Moss (Eds.), North Carolina State College, Raleigh, N.C.

Galerkin B.G., 1915, 'Series solution of some problems of elastic equilibrium of rods and plates (Russian) ', Vestn. Inzh.Tech., 19,897-908. Gasynski J. and G. Szefer, 1978, 'Axisymmetric problem of the punch for the

consolidating serni-space with mixed boundary penneability conditions', Archiwum Mechaniki Stosowanej, 30, 17-26. Ghaboussi I. and E.L. Wilson, 1973, 'Flow of compressible fluid in porous elastic media', International Journalfor Numerical Methodr in Engineering, 5,4 19-442.

Gibson R.E., K. Knight and P.W.Taylor, 1963, 'A critical expriment to examine theories of three dimensional consoüdation', Proceedings of European Conference on Soil Mechanics and Foundation Engineering, Wiesbaden, 1,69976.

Gibson R.E., R.L. Schifian, and S.L. Pu, 1970, 'Plane stralli and axially symmetric consolidation o f a clay layer on a smooth irnpervious base'. Quarterly Journal of Mechanics und Applied Mathematics, 23,505-520. Giraud A. and G. Rousset, 1996, 'Time-dependent behaviour o f deep clays', Engineering

Geology, 41, 181-195. Gladwell G.M.L.,1980, Contact Problems in the Theory of Elasticity, Sijthoff and Noordhoff,Alphen aan den Rijn, The Netherlands. Gladwell G.M.L.,1995, 'On contact problems for a medium with rigid flat inclusions of arbitrary shape', International Journal of Soli& and Structures, 32,383 -3 89.

Griffith A.A., 1921, 'The phenomena of rupture and flow in solids', Philosophical Transactions of Royal Socieîy of London, A22 1, 163- 1 97. Griffith A.A., 1924, 'The theory of rupture'. Proceedings of 1st International Congress in

Applied Mechonics, C. Bienzeno and LM.Burgen (Eds.), Wdtman, London, 55-63. Harding LW.and I.N. Sneddon, 1945, 'The elastic stress produced by the indentation of the plane surface o f a semi-infuite elastic solid by a rigid punch', Proceedings of Cambridge Philosophical Socieiy, 4 1, 16-26.

Harr M.E.,1966, Foundutiom of nieoretical Soi1 Mechanics, McGraw Hill, N.Y. Henshell R.D. and KG. Shaw, 1975, 'Crack tip fuiite elements are unnecessary', lnternational JourmI for Numerical Metho& in Engineering, 9, 495-507.

Hermann L.R.,1965, 'Elasticity equatioos for incompressible and nearly hcompressible materials by a variational theoremf, Journal of American Institute of Aeronuutics and Astronautics, 3, 1 896- 1900.

Hu J. and A.P.S Selvadurai, 1995, 'Influence o f tertiary creep on the uplift behaviour of a pipe embedded in a fiozen soil', Second International Confirence on Advances in

Underground Pipeline Engineering, Seattle, WA, I.K. Jeyapalan and M . Jeyapalan (Eds.), 345-358. Huang N.C. and S.G.Russell, 1985% 'Hydraulic fracturing of a saturated porous

medium-1:general theory', Theoret. Appl. Fracture Mech.. 4,20 1-2 1 3. Huang N.C. and S.G. Russell, 1985b, 'Hydraulic fiacturing of a saturated porous

medium-1:special cases', nieoret. Appl. Fracture Mech., 4,215222. Inglis C.E.,1913, 'Stresses in a plate due to the presence of cracks and sharp corners'. Transactions of the Institute of Naval Architects, 55,2 1 9-24 1 . Ingraffea A.R., 1 97%

Discrete Fracture Propagation in Rock: Laborutory Tests und

Finite Element AnaZysis, Ph.D. thesis, University of Colorado. U.S.A.

Ingraffea A.R., 1977b 'Nodal grafting for crack propagation studies', Internationol Journal for Numericol Methoh in Engineering, 1 1, 1 1 85- 1 1 87.

Ingraffea A.R. and T.J. Boone, 19 88, 'Simulation o f hydraulic fracture in poroelastic rock', Proceedings of Confireence on Numericol Methoh in Geomechanics, Swoboda (Ed.), Bakema,Rotterdam, Vol. 1,95- 105.

Invin, GR., 1957, 'Analysis of stresses and strains near the end of a crack travershg a plate', Transactions o f ASME,Journal afApplied Mechunics.

Invin, G.R., 1958, Fracture, in Handbuch der Physik, Vol. VI, Springer, Berlin.

Johnson K.L., 1 985, Contact Mechnics, Cambridge University Press, Cambridge, England.

Ju J.W., 1990, 'Isotropie and anisotmpic damage variables in continuum damage mechanics', Journal of Engineering Mechmics, ASCE, 116, No.12,276402769.

Kachanov L.M., 1958, ' T h e of ruphire process under creep conditions', Isv. A b d . Nuuk SSR Ofd. Teck ', 8,2603 1. Kassir M.K. and G.C. Sih, 1975, 'Three dimensional notch problems', in Mechanics of

Fracture, 2, G.C.Sih (Ed.), Noordhoff, Leyden, 173- 194. Kiyama T.,H. fita, Y. Ishijima, T.Yanagidani, K. Aoki, T. Sato, 1996, 'Permeability in anisotropic granite under hydrostatic compression and triaxial compression including p s t - failure region', Proceedings

of the 2nd North A merican Rock Mechanics Symposium:

NARMS' 96, Vol. 2, 1643-1650. Krajcinovic D. and GU. Fonseka, 1981, 'The continuous damage theory of brittle rnaterials', Journal of Applied Mechanics, 48,809-8 15. Krajcinovic D., 1984, 'Continuum damage mechanics', Applied Mechanics Reviews, 37, 1-6.

Kuroki T. Ito and K. Onishi, 1982, 'Boundary element method in Biot's linear

consolidation',Applied Mathematicul Mudelling, 6, 1 05- 1 10. Lambe T.W. and R.V. Whitman, 1969, &il Mechunics, John Wiley & Sons, New York.

Lan Q. and A.P.S. Selvadurai, 1996, 'Interacting indentors on a poroelastic half-space', Journal of Applied Mathematics and Physics, 47, 1 -22.

Lemaitre J. and J.L Chaboche, 1990, Mechanics of Solid .Uuteria&s,Cambridge University Press, Cambridge, England.

Lemaitre 1.. 1984, 'How to use damage mechanics', NucIem Engineering and Design, 80, 233-245.

Lewis R.W.and B.A. Schrefler, 1987, The Finite Elemenr .Clethod in the Deformaiion and Consolidation of Porous Media, John Wiley & Sons, Yew York.

Liebowitz, H. (Ed.), 1968, Fracture, Vol. 1-IV, Pergamoa Press. Oxford. Lur'e A I . , 1964, Three Dimensional Problems of the 7'heory of Elasticity, John Wiley, New York, N.Y.

Mandel J., 1950, 'Etude mathématique de la consolidation des sols. actes du colloque international de mécanique', Poitier, France, 4,9019.

Mandel J., 1953, 'Consolidation des sols (étude mathématique)', Geotechnique, 3, 287299.

McNarnee J. and R.E. Gibson, 1960% 'Displacement fimctions and linear transfomis appîied to diffusion through prous elastic media', Quart. J Mech. und Appiied Math.. 13,980 1 1 1.

McNamee J. and R.E. Gibson, 1960b, 'Plane strain and axiafly symmetric problems of the consolidation of a semi-infinite clay stratum', Quart. J Mech. and Applied Math., 13,

Melin S., 1991, 'On the directionai stability of wedging', International Journal of Fracture, 50,293-300. Muni V. and S . Valliappan, 1986, ' A universal optimum quarter point element',

Engineering Fracture Mechanics, 25(2), 237-258. Nay lor D.J., 1 974, ' Stresses in nearly incompressible materials by f ~ t elements e with

application to the caiculation of excess pore pressures'. Ihternationacll Journal for

Numerical Methoh in E w e e r i n g , 8,443-460. Ngo D. and A.C. Scordelis, 1967, 'Finite element analysis of reinforced concrete beams', Journal of Arnericun Concrete Inriif~ e 64, , 152- 163.

Nierode D.E., 1985, 'Cornparison of hydraulic fracture design methods to observed field

-

results',Journal of Petrulem Technology, 37, 183 1 1839.

Nilson A.H., 1968, 'Nonlinear analysis of reinforced concrete beams by the finite element method', AC1 Journal, 65, No. 9,757-766.

Owen D.R.J. and A.J. Fawkes, 1983, Engineering Fracture M e c h i c s : Numerical

Methoh and Applications, Pinendge Press Lirnited, Swansea, U.K. Palmer A.C. and J.R. Rice, 1973, 'The growth of slip surfaces in the progressive failure of overconsolidated clay', Proceedings of Royal Sociefy of London, Series A, 332, 527548,

Radok J.R.M.,1956, 'On the solution of problems of dynamic plane elasticity', Quarteriy of Applied Mathematics, 14,289-98.

Rashid Y .Ra,1968, 'Andysis of prestressed concrete pressure vessels', NucIear Engineering Design, 7,334-344. Rendulic L., 1936, 'Porenzîffer und Poren Wasserdruk in Tonen', Bauingenieur, 17, 559564.

Rice J.R. and M.P. Cleary, 1976, 'Some basic stress diffusion solutions for fluidsaturated elastic porous media with compressible constituents', Reviews of Geophysics and Space Physics, 14,227-24 1.

Rice J.R. and D.A. Simons, 1976, 'The stabilization of spreading shear faults by coupled

deformation-diffusion effects in fluid-infiltrated porous material', Journal of Geophysics Research, 81,53224334.

Rudnicki J.W., 1985, 'Effect of pore fluid diffusion on deformation and failure of rock', Mechanics of Geomaterials, 2. Bazant (Ed.), John Wiley & Sons, New York.

Ruina A., 1978, 'Influence of coupled deformation-diffusion effects on retardation of hydradic fracîure', Proceedings of 19th US.Symposium in Rock Mechnics, Nevada, USA, 274-282. Samaha H.R.and K.C.Hover, 1992, 'Muence of microcracking on the mass transport properties of concrete', ACI Materials Journal, Vol. 89, No.4,4 16424. Sandhu R.S. and E.L. Wilson, 1969, 'Finite element analysis of flow in saturated porous

media', Journal ofEngineering Mechanics Division, ASCE, 95 (EMJ), 641-652.

Sandhu RS., H. Liu and

K.J. Sin& 1977, 'Numerical performance of some

hite

element schemes for analysis of seepage in prous elastic media', Interncrtional Jourml of

Rock Mechmics and M ' h g Science and Geomechanics Abstrac~s,1, 177-194.

R., T. Stankowski, K.H.Gerstle and H.Y. Ko, 1983, 'Stress-shain curves for concrete under multiaxial load histones', NSF CME-80-0 1508, Department of Civil Scavuvo

Engineering, University of Colorado, Boulder. Schihan R.L., A.T-F Chen and J.C. Jordan, 1969, 'An analysis of consolidation theones', Journul of Soi1 Mechonies and Foundorion Divisions, ASCE, 95(SM 1 ), 285-

312.

Schifian R.L. and A.A. Fungaroli, 1965, 'Comolidation due to tangentid loads', Proceedings

4 6th

Inrernationd Conference on Soil and Foundulion Engineering,

Montreal, Canada, 1 , 1 88- 192. Schrefler B.A. and L. Sirnoni, 1987, 'Non-isothermal consolidation of unbounded porous media ushg mapped infinite elements', Comm. Appl. Numer. Meth., 3,445-452.

Selvadurai A.P.S., 1979, EZmtic Anabsis of &il-Foundation hteraction; Developments in Geoteehnical Engineering, Vol. 17, Elsvier Scientific Publications, The Neiherlands.

Selvadurai A.P.S. and B.M. Singh, 1984, 'On the expansion of a penny-shaped crack by a rigid circdar dix: inclusion', Intematio~i&Nimd ojFmcrwe, 25,69-77. Selvadurai A.P.S., 1985, 'On an integral equation goveming an intemally indented penny-shaped crack', Mechortics Resewch Communications, 12,347-35 1.

Selvadurai A.P.S. and K.R. Gopal, 1986, 'Consolidation analysis of the screw plate test', Proceedings of 39th Canadian Geotechnical Society Conference, Ottawa, Canada, 167178. Selvadurai A.P.S. and R. Karpurapu, 1989, 'A composite infinite element for modelling of unbounded saturated soi1 media', Journal of Geotechnical Engineering, ASCE, 115, 1633-1646. Selvadurai A.P.S. and M.C. Au, 1989, 'Crack behaviour associated with contact problems with non-linear interface constraints', in Boundury element techniques: Applications in Engineering, C.A. Brebbia and N.G. Zamani (Eds.), Computational Mechanics Publications, Boston, 1- 17. Selvadurai A.P.S. and M.C. Au, 1991, 'Damage and visco-plasticity effects in the indentation of a polycrystalline solid', Proceedings of Plasticity 91: 3rd International Symposium on Plasticity and Its Current Applicatiom, J.P.Boehler and A.S. Khan (Eds.), Elsevier Applied Science, London, 405-408. Selvadurai A.P.S., 1994a, 'A unilateral contact problem for a rigid disc inclusion embedded between two dissimilar elastic ha1f-spaces', Quarterly Journal of Mechanics and Applied Mathematics, 47,493-5 10. Selvadurai A.P.S., 1994b, 'Separation at a pre-fractured bi-material geological interface', Mechanics Research Communications, 2 1,83088. Selvadurai A.P.S., 1994c, 'Transverse loading of unidirectionally reinforced composites: the role of maûix damage and matrix fracture', in Dumage Mechanics of Composites, ASME 1994, D.H. Allen and J.W. Ju (Eds.), AMD-Vol. 185, 117-127.

Selvadurai A.P.S. and Z.Q.Yue, 1994, 'On the indentation of a poroelastic layerl, International Journal for Numericol and Analytical Methoh in Geomechanics, 18, 161175. Selvadurai A.P.S. and T.S. Nguyen, 1995, 'Computational modelling of isothermal consolidation of hctured porous media', Computers and Geotechnics, 17, 39-73. Selvadurai A.P.S. and A. ten Busschen, 1995, 'Mechanics of the segmentation of an embedded fiber, Part II: Computational modelling and comparisons', Journal of Applied Mechanics, Transactions ofASME. 62,98- 107.

Selvadwai A.P.S. (Ed.), 1996, Mechanics of Poroelastic Media, Kluwer Academic Publ., A. A. Dordrecht, The Netherlands.

Selvadurai A.P.S.,

1997% 'Certain non-classical effects in contact mechanicsl,

Proceedings of Contact Mechanics ' 97, Madrid, Spain (in press). Selvadurai A.P.S., 1997b, 'The elastic cornpliance of a punch bonded to a weakened elastic half-space', (in preparation). Sidoroff F., 1980, 'Description of anisotropic damage application to elasticity', Int. The. Appl. Mech. Symposium Senlis, France, Physical Non1inearities in Structural A nalysis. J. Hult and J. Lemaitre (Eds.), Springer-Verlag, Berlin Heidelberg, New York, 237-244.

Shephard M.S., N.A.B. Yehia, G.S. Burd and T.J.Weidner, 1985, 'Automatic crack propagation tracking', Computers and Structures, 20(1), 2 11-223. Shih C.F.,H.G. de Lorenzi and M.O.German, 1976, 'Crack extension modelling with singular quadratic isoparametric elements', International Journal of Fracture, 12, 647651.

Shih C.F., H.G.de Lorenzi and H.R. Andrews, 1979, 'Studies on crack initiation and stable crack growth', in Elostic-Plastic Fracturing, LD. Landes, J.A. Begley and G.A. Clarke (Eds.), ASTM STP 668,65- 120. Shiping L., L. Yushou., L. Yi, W. Zhenye, and 2. Gang, 1994, 'Permeability-Sûain epuations corresponding to the complete stress-strain path of Yirizhuang Sandstone',

International Journai of Rock Mechanics and Mining Science and Geomechanics Abstracts, 31,383-391. Sih G.C., 1974, 'Shain energy density factor applied to mixed mode crack problems',

Inter~fional Journal of Fracture, 10,3 05-32 1. Sih G.C.,1991, Mechanics of Fmcture Initiation and Propagation, Kiuwer Academic

Publisher, Dordrech~The Netherlands.

Sih G.C. and Theocaris, P.S. (Eds.), 1979, MLred Mode Crack Propagation, Proceedings of 1st US.-Greece S'mposiurn, Athens, Greece, Sijthoff and Noordhoff Alpen aan den Rijn, The Netherlands. Simo

J.C.and J. W. Ju, 1987, 'Strain- and stress-based continuum damage models-II.

-

Computational aspects', Internatonol Journal of &Iids und Structures, 23, No. 7, 84 1

869.

Sirnous D.A.,1977, 'Boundary layer andysis of propagating mode II crack in porous

elastic media', Journai of Mechanics and Physics ofSoi&, i5,99-115.

Skempton A.W., 1954, 'The pore pressure coefficients A and B', Geotechnique, 4, 143147.

Smith M.B.,1985, 'Stimulation design for short, precise hydraulic fractures', Society 4 Petroleum Engineers Journal, 25,37 1-3 79. Smith LM. and D.V. Grifiths, 1988, Progrunrming thefinite element method, John Wiley & Sons, New York.

Smith R.N.L. and J.C. Mason, 1982, 'A boundary element rnethod for c w e d crack

problems in two dimensions', Boundary Element Methoris in Engineering, C.A. Brebbia

(Ed.), Springer-Verlag, Berlin, 472-484. Sneddon LN., 1946,' The distribution of stress in the neighborhood of a crack in an efastic solid', Proceedings of the Royal Socie& Series A, 187, 220-260. Spooner D.C. and J.W. Dougill, 1975, 'A quantitative assessrnent of damage sustained in

concrete during compressive loading', Magazine of Concrete Research, 27, 1 5 1- 160.

Szefer G. and J. Gaszynski, 1975, 'Axisymmetric punch problem under condition of consolidation',Archiwum Mechaniki Stosowanej, Warsaw, Poland, 27,497-5 15. ten Busschen A. and A.P.S. Selvadurai, 1995, 'Mechanics of the segmentation of an

embedded fiber, Part I: Experimental investigations', J o d of Applied Mechmies, Trmactions of ASME, 62, 87-97.

Terzaghi K., 1923, 'Die Berechnung der durchlassigkeitsziffer des tones aus dem verlauf der hydrodynamischen spannungserscheinungen', Ak. der Wissenschafien in Wien, Sitzungsberichte mathematisch-naturwissenschaftliche Klasse. part 114 1 32(3/4), 125138-

Tong P. and T.H.H. Pian, 1973, 'Onthe convergence of the f ~ t element e method for problems with singulariîy', International Journal of Solith und Structures, 9,3 13-32 1.

Valliappan S. and V. Murti, 1985, 'Automatic remesbg technique in quasi-static and dynamic crack propagation', Proceedings of NUMETA '85 Confieence, Swansea, U .K., 107-116. Vandsmme L., E. Detoumay and A.H.D. Cheng, 1989, 'A WO-dimensional poroelastic

displacement discontinuity method for hydraulic fracture simulation', International Journal for Numerical and Analytical Methoh in Geomechnics, 13,2 15-224. Ve-t

A., 1965, Discussion, SUth hterncrtional Conference on Soil Mechanics and

Foundution Engineering, 3,40 1-402.

Wawrzynek P.A. and A.R. Ingraffea, 1989, 'An interactive approach to local remeshing

around a propagating crack', Finite Eiemenf Amlysis and Design, 5, 87-96. Westergaard H.M.,1938, 'Bearing pressures and cr;ickst,Journal of Applied Mechanics. 6,49053.

Williams M.L., 1957, 'Onthe stress distribution at the base of a stationary crack'. Journal of Appkd Mechanics, ASME, 24, 10%1 14. Wolfram Research Inc., 1993, MATHEMATICA, Version 2.2, Champaign, Illinois, USA.

Yoffe E.H.,195 1, 'The moving Griffith crack', Philosophical Magazine, 42, 739-750.

Yue Z.Q.and A.P.S. Selvadurai, 1994, 'On the asyrnmeûic indentation of a consolidating poroelastic half space', Applied Mathematical Modelling, 18 (4), 170- 185.

Yue Z.Q.and A.P.S. Selvadurai, 1995, 'Contact problem for saturated proelastic solid', Journal of Engineering Mechunks, ASCE, 121 (4), 502-5 12.

Zaretskii Y.K., 1972, 'Theory of Soil Consolidation, (Translated fiom Russian) ', U.S.

Department of Commerce National Technical Information Service, Verginia, TT7050131.

Zienkiewicz O C ,1977, The Finite Eiement Method, McGraw Hill, New York. Zoback M.D.and J.D.Byerlee, 1975, 'The effect of microcrack dilatancy on the permeability of Westerly Granite',Journal ojGeophysica1 Reseurch, 80,752-75 5.

FINITE ELEMENT FORMULATIONS

OF POROELASTIC MEDIA

The generai formulation of finite element analysis of ihreedhensional consolidation is derived without refemng to any specific types of element. The goveming equations considering the body forces take the form:

A boundary value problem has to be solved which requires that the goveming equations

B. 1 and B.2 are satisfied within the domain R and the appropriate boundary conditions are satisfied on the boundary B of the domain. Using Galerkin's technique, the governing

equations can be transfonned into matrix equations w h m unknowns are the nodal

displacements and pore pressures. It has k e n shown by Sandhu and Wilson (1969). and

Ghaboussi and Wilson (1973) that to ensure stability of the solution, the nodal displacements are assigned an order different mgher) to the stresses and pore pressures. Letting u, (i=1,2,3; hl,...,N ) be the nodal displacements for N discrete points in R and pK (K=l,

..., n) the nodal pore pressures in n nodes

at

an arbitrary time

t.

The

displacement vector and pore pressure for any arbitrary point with coordinates x.J in the dornain R are approximated by the following dations:

Where N: and N i are, respectively, shape functions for the displacement field and for the

pore pressure field; &1, ..., N and K=I , ..., n; u, is the displacement of the soi1 skeleton at node J in the ith direction. The indices

in capital letten (e.g. J and K) refer to nodal

values, while the indices in small letters (e.g. i andj) refer to coordinate directions. Also summation convention is adopted. In general N: and N [ can be different but both and

Ni

B.1

Galerkin Formulation for the Equilibrium Equatioa

Si

must exhibit CO continuity.

Applying the Galerkin's weighted residual method to the equilibrium equation B. 1 results

in following

Where I= 1, ...,N. Using Green's theorern, Equation B.5 becornes

Where B is the boundary of R. Substituthg the approximation equations B.3 and B.4 into the above equation, one obtains

Where fiom the constitutive equations 2. la and strain tensor equation 2.2

Equation B.7 can be written in matrix forrn as

Where {6} and ( p } are the vectors of nodal displacements and pore pressures:

(B.10)

(F,) and {Ft} are respectively the vectors representing the body force and the traction applied at the boundary B. The components of the coupling matrix [Clwhich results from the interaction between the soi1 skeleton and the pore fluid are given by the following equation

The stiffness matrix [KIof the soil skeleton is generally written in the following form:

(B.12)

where [Dlis the stress-strain rnatrix for the soil skeleton. For an isotropic linear elastic material [Dl depends on two elastic constants p and

)c

(or E and v ) . [BI is the matrix

relating strains to nodal displacements which depends on the shape iùnctions N ; .

B.2 Galerkin Formulation for the Fluid Continuity Equation Using concept of Gaierkin method in association with flow continuity equation B.2, we

obtain the weighted residual equivalent of the Equation B.5 as

(B.13)

Applying Green's theorem to Equation B. 13, one obtains

(B.14)

Substituting B.3 and B.4 into the above equation, one obtains

The Equation B. 15 c m be written in matrix fonn as

where the permeability matrix [Hltakes the form:

(B.17)

and compressibility rnatrix [El of pore fluid takes the fom:

(B.1 8)

and {F } is the inward fluid flux through the boundary B. The matrix equations B.9 and q

B. 16 are the fuiite element approximations for the goveming equations of poroelasticity. 8.3 Time intcgntioo and Stabiüty

In order to mode1 the nonlinear behaviour of materials such as damage or plasticity phenornena, it is necessary to use the governing equations in an incremental form (Zienkiewicz et ai., 1977; Lewis and Schrefler, 1987). By differentiating the equiiibrium equation B.9, the system of incrementai coupled equations is obtained as follows

(B.19)

The matrices in above equation relate increments of responses (a, dp) to increments of extemal driving forces (e.g. dF3. These matrices are dependent of the curent state (6, p)

of the system. The system of coupled equations B.19 is discretized in the tirne domain by foilowing

finite difference scheme for any variable X in the system:

O

1

where At is time increment, X ,X ,T are the values of quantity Xat different times t, t+ At and r + y b respectively; y is a value between O and 1. When ~ 0the, finite diflerence

scheme is called hlly explicit; when y =1, it is called fully implicit; when y =OS, it is cdled the Crank-Nicholson scheme. Applying this finite difference scheme to &/dt, and dpldf in Equations B. 19, results in the

following elemental matrix equation:

where

F = force vectors due to extemal tractions, body forces and flow field; u, p, = nodal displacements and pore pressures at time t;

& = t h e Uicrement and y is the time integration constant.

M e n u, and p, are known, the solution of Equation B.22 results in ut+, and pl+,,. The

time integration constant y varies between O and 1. The criteria for the stability of the integration scherne given by Booker and Small (1976) require that y 2112. According to Lewis and Schrefler (1 987) and Selvadurai and Nguyen (1995), the stability of solution

can generally be achieved by selecting values of y close to unity.