Large strain primary and secondary consolidation of a soft clay using the finite element method

MISSING IMAGE

Material Information

Title:
Large strain primary and secondary consolidation of a soft clay using the finite element method
Physical Description:
xiv, 191 leaves : ill. ; 29 cm.
Language:
English
Creator:
Ahmed, Zafar, 1967-
Publication Date:

Subjects

Subjects / Keywords:
Soil consolidation -- Mathematical models   ( lcsh )
Finite element method   ( lcsh )
Soil mechanics -- Mathematical models   ( lcsh )
Differential equations, Nonlinear   ( lcsh )
Civil Engineering thesis, Ph.D   ( lcsh )
Dissertations, Academic -- Civil Engineering -- UF   ( lcsh )
Genre:
bibliography   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph.D.)--University of Florida, 1999.
Bibliography:
Includes bibliographical references (leaves 183-190).
Statement of Responsibility:
by Zafar Ahmed.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 030470892
oclc - 42680201
System ID:
AA00013546:00001


This item is only available as the following downloads:


Full Text










LARGE STRAIN PRIMARY AND SECONDARY CONSOLIDATION
OF SOFT CLAY USING THE FINITE ELEMENT METHOD
















By

ZAFAR AHMED


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY


UNIVERSITY OF FLORIDA

1999


* ..... : :: *: ..


iiii ;..:ll iiiE. :i:
:i: .:iii; ;i :: i::. .

iiiiii~i~i iiEii~ii!i~ii ii.:
.,iiiiiii;; i~iii!i: ': i%






























To my parents and Runu










ACKNOWLEDGMENTS

First, I would like to thank Dr. M. C. McVay for his overall assistance and

guidance during the completion of this research project. His deep insight and intuition in

the vast arena ofgeotechnical research always made a seemingly impossible task worth

pursuing. My deep reverence goes for the endeavor, spirit, and countless hours he

dispensed to make things work, to straighten out any bottleneck.

I also express gratitude to the members of my committee Dr. F. C. Townsend,

Dr. J. L. Davidson, Dr. D. G. Bloomquist and Dr. W. W. Hager for their guidance and

suggestions.

I am indebted to all the friends I made during the course of my study at the

University of Florida. Special thanks go to Kailash, Paulo and Pedro for their moral

support, help, and most of all their invaluable camaraderie. I had many happy hours with

them.

Finally, I extend gratitude to my parents, who were always an impetus for my

education.










TABLE OF CONTENTS

age

ACKNOWLEDGMENTS ....................................................... .............. iii

LIST OF TABLES ...................................................................... .......... iii

LIST OF FIGURES ............................................................................. x

A B STR A C T ......................................... ................................. ......... ... xiii

CHAPTERS

1 INTRODUCTION ..........................................................................

1.1 Introduction .................................. .....................................1
1.2 Research Focus .................................................................... 2
1.3 Scope of Work ..................................... ..................................5

2 LITERATURE REVIEW ................................................................ 7

2.1 Consolidation ...................................................................... 7
2.2 Large Strain Multiplicative Plasticity ............................................. 9

3 THEORY OF MIXTURES ............................................................... 12

3.1 Basic Theory ........................................................................ 12
3.2 K inem atics .......................................................................... 12
3.3 Average Quantities ................................................................. 15
3.4 Balance Laws ..................................................................... 17
3.4.1 Balance of Mass ........................................................... 18
3.4.2 Balance of Linear Momentum .............................................. 19
3.4.3 Balance of Angular Momentum ............................................20
3.4.4 Balance of Energy ............................................................21
3.4.5 Entropy Inequality ......................................................... 23
3.5 Balance Laws for Saturated Soil ................................................. 25
3.5.1 Balance of Mass ............................................................ 25
3.5.2 Balance of Linear Momentum ...........................................25
3.5.3 Balance of Energy ........................................................... 29






3.5.4 Entropy Inequality ......................................................... 32

4 VARIATIONAL EQUATIONS, CONSTITUTIVE THEORIES ..................34

4.1 Boundary Value Problem.......................................................... 34
4.2 Strong Form ....................................................................... 35
4.3 Weak Form .............................................................. 37
4.4 Time Descretization Scheme ...................................................... 39
4.5 Large Deformation Plasticity Model for Soil Skeleton ...................... 42
4.6 Constitutive Law for Fluid Flow ............................................. 49

5 HYPERELASTIC-PLASTIC MCC MODEL ............................................ 51

5.1 Introduction .................................................... .... ... .......... .. 51
5.2 Hyperelastic Model ...................................... ...................... .... 52
5.3 Plasticity Model .................................. ................................. 56
5.4 Hardening Law.................................................................... 57
5.5 Flow Rule ........................................... ......................61
5.6 Consistent Tangent Moduli ...................................................... 63

6 HYPERELASTIC-VISCOPLASTIC MCC MODEL ................................... 69

6.1 Introduction ......... .. ................. ........................................ 69
6.2 Flow Rule ............ ................ ............................... ..... .. 70
6.3 Consistent Tangent Moduli ...................................................... 73

7 LINEARIZATION ............................................ ........................... 76

7.1 Prelim inaries ...................... ........................................... ..... .. 76
7.2 Linearization of Strong Form ................................................... 77
7.2.1 Equation of Equilibrium................................................. 78
7.2.2 Equation of Flow Continuity ............................................... 81
7.3 Linearization of Weak Form ...................................................... 85

8 FINITE ELEMENT FORMULATION ............................................... 88

8.1 Finite Element Framework ....................................................... 88
8.2 Matrix Equations ................................................................. 89

9 NUMERICAL EXAMPLES ............................................................. 97

9.1 Mixed Element ..................................................................... 97
9.2 One-dimensional Hyperelastic Consolidation........ ............................. 98
9.3 Plane Strain Hyperelastic Consolidation ...................................... 104







10 POLK COUNTY EXPRESSWAY ..................................................... 109

10.1 Phosphatic Waste Clay ......................................................... 110
10.2 Wick Drain and Instrumentation ................................................ 116
10.3 Constitutive Model Parameters .................................................. 118
10.4 FE Mesh ........................................................................... 129
10.5 Prediction of Primary Consolidation ......................................... 132
10.6 Prediction of Secondary Consolidation ..................................... 140

11 CONCLUSIONS .......................................................................... 145

APPENDICES

A MATHEMATICAL DERIVATIONS ................................................... 149

A. 1 Gradient of the Jacobian, J...................................................... 150
A.2 Balance of Energy of Saturated Soil .............................. .......... 150
A.3 Weak Form of DIVP+pg = 0 .............................................. 151
A.4 Weak From of div v + div = 0 .............................................. 152
A.5 Area Transformation of Flow Rate ........................................... 153
A.6 Additive Decomposition of Principal Natural Strain .............. ........ 154
A.7 Proof of Piola Identity: DIV Y = J div y .................................. 155
A.8 Linearization ofF, F ........................................................... 156
A.9 Linearization of J, j ........................................................... 157
A.10 Linearization of po ..............................................................158
A.11 Relation between Tensors A and D ......................................... 158
A. 12 Relation between Tensors a and d .......................................... 159
A. 13 Spectral Decomposition of b, C .................................................. 159
A. 14 Derivation of 6eA/ .C .............................................................. 162
A. 15 Derivation of 8aM1^/C ......................................................... 162
A.16 Push Forward of aM A)/OC ..................................................... 164
A.17 Variation of k, K ............................................................... 166
A.18 Variation of grad 0 ................................................................ 167
A.19 Variation of ............................................................ 167
A.20 Variation of grad 'q: ............................................................. 168
A.21 Variation of GRAD q: P ........................................................ 170
A.22 Variation of grad\y ........................................................... 171
A.23 Variation of GRAD ~y V ....................................................... 173
A.24 Hand Calculation of One-dimensional Large Strain,
Hyperelastic Consolidation ................................................ 174







B FINITE ELEMENT MATRICES ...................................................... 176

C LABORATORY CONSOLIDATION TEST DATA ................................ 180

REFERENCES ...................................... ......... .................. ......... ..... 183

BIOGRAPHICAL SKETCH .................................................................... 191










LIST OF TABLES


Table page

4.1 Genearlized trapezoidal method .............................................41

10.1 Summary of laboratory consolidation test results
(initial exploration) ..................................................... 114

10.2 Summary of laboratory consolidation test results
(later exploration) ....................................................... 114

10.3 Hyperelastic-plastic MCC model parameters for large
strain simulation of laboratory consolidation tests ................ 125

10.4 Hyperelastic-plastic MCC model parameters for small
strain simulation of laboratory consolidation tests .............. 125

10.5 Hyperelastic-plastic MCC model parameters for
simulation of laboratory consolidation test ..................... 126

10.6 Hyperelastic-viscoplastic MCC model parameters for large
strain simulation of laboratory consolidation tests ................. 128

10.7 Settlement cell/plates represented by 2.44 m deep pond .................130

10.8 Settlement cell/plates represented by 4.57 m deep pond ................ 130


10.9 Settlement cell/plates represented by 7.62 m deep pond .................. 130

10.10 Material parameters for hyperelastic-plastic consolidation
(pond depth 2.44 m) .................................................. 133

10.11 Material parameters for hyperelastic-plastic consolidation
(pond depth 4.57 m) ................................................. 134

10.12 Material parameters for hyperelastic-plastic consolidation
(pond depth 7.62 m) .................................................. 134






10.13 Comparison of primary consolidation settlements at 485 days:
large strain versus small strain ........................................137

10.14 Location ofpeizometers ..................................................... 138

10.15 Material parameters for hyperelastic-viscoplastic consolidation
(pond depth 2.44 m) ................................................ 141

10.16 Material parameters for hyperelastic-viscoplastic consolidation
(pond depth 4.57 m) ................................................ 141

10.17 Material parameters for hyperelastic-viscoplastic consolidation
(pond depth 7.62 m) ................................................ 142

10.18 Comparison of creep settlements of different ponds ..................... 144











LIST OF FIGURES


Figure Page

3.1 Geometric representation of kinamatics
for a two-phase mixture ................................................. 13

3.2 Balance of mass: flow through a control volume ......................... 18

3.3 Balance of mass: total mass of material point X
is not conserved in 4,(X) ............................................... 28

4.1 Prescribed boundary conditions of spatial domain ..................... 35

4.2 Illustration of multiplicative decomposition of the
deformation gradient ................................................... 43

5.1 Yield surface of the MCC Model in p-q plane .............................. 56

5.2a Unilogarithmic compressibility law ............................................ 58

5.2b Bilogarithmic compressibilty law............................................ 58

5.3 Limit of validity: comparison between unilogarithmic
and bilogarithmic compressibility laws ................................ 59

9.1 D9P4 mixed element ................................. ......................... 98

9.2 FE mesh and initial pore water pressures for
one-dimensional consolidation problem .............................. 99

9.3 One-dimensional hyperelastic consolidation: variation
of total potential with time ........................................... 100

9.4 One-dimensional hyperelastic consolidation: isochrones
of constant Cauchy pore pressure ................................. 102

9.5 One-dimensional hyperelastic consolidation: variation
of average degree of consolidation with time factor ..................... 103







9.6 FE mesh for plane strain hyperelastic consolidation example .............105

9.7 Plane strain hyperelastic consolidation: variation of centerline
excess pore pressure at depth z = a with time .............. ... 106

9.8 Plane strain hyperelastic consolidation: isochrones of
constant Cauchy pore pressure .................................... 108

10.1 SPT boring logs for tests reported in Table 10.1 (Source: PSI) ......... 112

10.2 SPT boring logs for tests reported in Table 10.2 (Source: PSI) .......... 113

10.3 Permeability versus void ratio plot from CRS
consolidation test ............................. ...................... .. 115

10.4 Void ratio versus effective plot from CRS
consolidation test ............................. ...................... 116

10.5 Plan of wick drain installation ................................................ 117

10.6 Principle of consolidation with wick drains .............................. 117

10.7 Instrumentation plan: surcharge area no. 1 of Polk County Expressway,
Section 5 (Source: Atlanta Testing & Engineeing) ............... 119

10.8 Instrumentation plan: surcharge area no. 2 of Polk County Expressway,
Section 5 (Source: Atlanta Testing & Engineeing) .............. 120

10.9 Cross-section of clay slime ponds of surcharge area no. 1
(Source: PSI) .................... ............ ........... ... ...... .... 121

10.10 Cross-section of clay slime ponds of surcharge area no. 2
(Source: PSI) ............................... ............ ......... ... .. 122

10.11 FE mesh for oedometer cell ................................................... 123

10.12 Large strain, hyperelastic-plastic simulation
of laboratory consolidation tests ................................. 123

10.13 Small strain, hyperelastic-plastic simulation
of laboratory consolidation tests .................................. 124

10.14 Hyperelastic-plastic simulation of laboratory






consolidation test: large strain versus small strain ............... 126

10.15 Large strain, hyperelastic-viscoplastic simulation
of laboratory consolidation tests .................................... 127

10.16 Schematic of contributive cylinder of soil
surrounding wick drains ............................................ 131

10.17 FE meshes for different pond depths ......................................... 132

10.18 Hyperelastic-plastic consolidation settlement:
pond depth 2.44 m ................................................... 135

10.19 Hyperelastic-plastic consolidation settlement:
pond depth 4.57 m ................................................... 136

10.20 Hyperelastic-plastic consolidation settlement:
pond depth 7.62 m ................................................... 136

10.21 Piezometer data versus prediction of total Cauchy
pore pressure: pond depth 4.57 m ................................... 138

10.22 Piezometer data versus prediction of total Cauchy
pore pressure: pond depth 7.62 m ................................ 139

10.23 Hyperelastic-viscoplastic consolidation settlement:
pond depth 2.44 m .............................. .................... 142

10.24 Hyperelastic-viscoplastic consolidation settlement:
pond depth 4.57 m .................................................. 143

10.25 Hyperelastic-viscoplastic consolidation settlement:
pond depth 7.62 m .................................................. 143









Abstract of Dissertation Presented to the Graduate School of the
University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

LARGE STRAIN PRIMARY AND SECONDARY CONSOLIDATION
OF SOFT CLAY USING THE FINITE ELEMENT METHOD

By

Zafar Ahmed

August 1999

Chairman: Michael C. McVay
Major Department: Civil Engineering

A mathematical framework for large strain consolidation of fully saturated

soil is presented in this study. Strong and weak forms of the boundary-value problem

are derived using both the material and spatial descriptions. The algorithmic treatment

of large strain of the solid phase is based on multiplicative decomposition of the

deformation gradient and is coupled with the algorithm of fluid flow via the

Kirchhoff pore water pressure. Balance laws for saturated soil are written following

the classical theory of mixture considering the soil-water mixture as a two-phase

continuum. Since the motion of the fluid phase only affects the Jacobian of the solid

phase, balance laws are derived completely in terms of the motion of the soil

skeleton.

Both geometric and material nonlinearities are incorporated in the stated

mathematical model. Geometric nonlinearity is represented by the Jacobian and its

variation along with proper stress and strain characterization based on the updated


xiii






configuration. Material nonlinearities of soil skeleton are represented in this study by

two different constitutive relations (elastoplastic and viscoplastic) based on a

modified Cam-Clay (MCC) model of critical state soil mechanics.

The proposed mathematical model is linearized in order to obtain consistent

tangent operators and subsequently implemented in the displacement-based finite

element code PlasFEM, developed by the Geotechnical Engineering group at the

University of Florida. Numerical predictions of primary (consolidation and swell) and

secondary (creep) consolidation settlement are made using PlasFEM for phosphatic

waste clay, found at the construction site of the Polk County Expressway in Lakeland,

Florida, in the presence of vertical wick drains and subject to surcharge loading,

unloading, and subsequent reloading (i.e., road construction). Simulation results

showed a very good agreement with the field measurements.


xiv












CHAPTER 1
INTRODUCTION



1.1 Introduction

Non-linear response ofgeotechnical structures typically results from plastic

yielding and large deformation of the soil skeleton. There are many classical geotechnical

applications where non-linear effects due to these two factors could critically influence

the outcome of a numerical analysis. Two examples are large movement of slopes and

tilting of a tower due to P-8 effect. The impact of large deformation and plastic response

is most evident in soft clays. It is a well-known characteristic of clays that considerable

time is required for the occurrence of the compression caused by a given increment of

load. Two phenomena contribute to this large time lag. The first is due to time required

for the escape of the pore water. It is called the hydrodynamic lag or consolidation, a

phenomenon which involves transient interaction between the solid and fluid phases of a

soil-water mixture. The second phenomenon is called plastic time lag or secondary

compression. The slow continued compression that continues after the excess pore

pressures have substantially dissipated is called secondary compression. Secondary

compression occurs because the relationship between void ratio and effective stress is

usually somewhat time-dependent: the longer the clay remains under a constant effective

stress, the denser it becomes.








Various problems of coupled fluid flow and deformation in porous media arise

frequently in the fields of geotechnical engineering, groundwater hydrology and plate

tectonics, among others. Slope stability analysis, surcharge loading for consolidation

drainage, embankment, excavation, settlement of bridge approach pavements etc. are a

few examples of a wide spectrum of geotechnical applications that involves

consolidation.

The behavior of soil during consolidation is governed by the differences between

the total stresses acting on the soil mass and the pore pressure. In most practical field

cases it is necessary to describe the effective stressfield to characterize the strength and

deformation properties of the soil. Although the no-flow (undrained) and the free-flow

(drained) conditions can be analyzed using a single-phase continuum formulation,

consideration of a two-phase soil-water relationship in a saturated soil medium is

essential in characterizing the soil behavior during the transient period of excess pore

pressure dissipation. The physics involved in consolidation phenomenon requires that

appropriate numerical analysis address coupled response of solid and fluid phases.



1.2 Research Focus

The mathematical structure and numerical analysis of nonlinear consolidation at

small strains are fairly well developed and adequately documented [1-10]. The general

approach is to write the linear momentum and mass balance equations in terms of the

solid displacement and fluid potential (or pore water pressure), and then solve them

simultaneously via a two-fold mixed formulation. The small strain assumption simplifies

the linear momentum balance equation since it produces an additive form of elastic and








plastic deformations. In the context of finite element analysis, the small strain assumption

also simplifies the mass conservation equation since the volume change of the mixture

becomes a linear function of the nodal solid displacements.

In spite of substantial development of computational methods for small strain

consolidation, mathematical models capable of handling the problem of coupled fluid

flow and large deformation of the soil matrix are not developed well enough to be useful

for routine analysis of prototype geotechnical structures. Extensions of the small strain

formulation of the classical consolidation equations to large deformation are based

primarily on the use of rate-constitutive equations [8,9,12-14]. In addition to the

restrictions of small elastic strains imposed by this hypoelastic formulation, it also

obscures a proper definition of'mean gradient' and 'average volume changes' necessary

for imposing the mass conservation equation at finite increments. Consequently, second-

order terms in the hypoelastic extension are ignored, particularly in the mass conservation

equation, which leads to a degradation of accuracy when the load increment is large.

The present study adopts an alternative formulation for large strain elastoplasticity

based on the multiplicative decomposition of the deformation gradient. This method

completely circumvents the 'rate issue' in large deformation analysis [17,18], and allows

for the development of large elastic strains. Multiplicative decomposition technique

better represent the particulate nature of soil, much like for metals from its crystal

microstructure. It provides a means for describing mathematically the relationships

between the reference configuration, the current configuration, and the unloaded, stress-

free intermediate configuration of a soil assembly subjected to large deformation in the

microscopic sense. A more recent development [19,20] indicates that the multiplicative






decomposition technique can be exploited to such an extent that the resulting algorithms

may inherit all the features of the classical model of small strain plasticity.

Proper characterization of fluid flow is another long-standing issue in large

deformation consolidation analysis. Classical theory of mixtures [21-26] is employed in

this study to describe coupled response of solid and fluid phases. Accordingly soil-water

mixture is viewed as a two-phase continuum, appropriate balance principles that govern

the interaction between the solid and fluid constituents are derived. In contrast to

previous formulations of the mixture theories, however, this study focused on the motion

of the solid phase alone and uses the constitutive flow theory in terms of relative motion

of the fluid with that of the solid [27]. Spatial form of generalized Darcy's law is used to

describe the constitutive relation of the fluid phase.

Due to its simplicity and practicality, modified Cam-Clay (MCC) model [28, 29]

of critical state soil mechanics is adopted in this study to represent the nonlinear

responses (plastic and viscoplastic) of the solid phase. Important features like

hyperelasticity and bilogarithmic compressibility law are added to MCC model to better

simulate the behavior of soft clay in large strain regime.

Finally the mathematical framework of non-linear large deformation

consolidation analysis is implemented in a displacement-based finite element code

PlasFEM, which is been developed by the Geotechnical Engineering group at the

University of Florida over last few years. Numerical analysis is performed to study the

consolidation behavior of a field case of low solid-content phosphatic waste clay.






1.3 Scope of Work

The dissertation is organized in eleven chapters. Chapter 2 presents literature

review and historical development of the theory of mixtures for porous media based on

balance laws, micromechanical approach for large strain based on the theory of

multiplicative decomposition. In Chapter 3, general kinematics and balance laws for a

non-interacting mixture of non-polar constituents are derived in the light of modern

theory of mixtures. Balance equations, specific for a fully saturated soil media (two-

phase), are further reduced from generalized equations. In Chapter 4, field equations or

strong form of the boundary-value problem of consolidation phenomenon are established

from balance equations, corresponding weak forms are derived for use in subsequent

finite element formulation. Constitutive theories for solid and fluid phases, appropriate

for large strain, are outlined. Constitutive models for phosphatic waste clay (low solid-

content clay) are discussed in Chapters 5 and 6. Hyperelastic-plastic MCC (modified

Cam-Clay) model, suitable for primary consolidation response, is presented in Chapter 5.

Chapter 6 presents the development of a hyperelastic-viscoplastic MCC model, used for

simulation of time-dependent secondary compression. Explicit expressions for the

consistent tangent moduli of the constitutive models are derived in the framework of

large deformation theory, based on multiplicative decomposition as mentioned before. In

Chapter 7, corresponding variational equations for boundary value problems are

developed and linearized for implementation into a finite element code. Chapter 8

presents the implementation issue of the governing equations, i.e., matrix formulation for

finite element code. Chapter 9 presents numerical examples one and two-dimensional

hyperelastic consolidation. These examples demonstrate significance of large strain on









consolidation responses compared to the same for small strain formulation. Chapter 10

presents a study of consolidation phenomena of phosphatic waste clay, deposited at the

construction site of Polk County Expressway (a multi-lane expressway around Lakeland,

Florida). Numerical simulations are run for cases of hyperelastic-plastic (primary) and

hyperelastic-viscoplastic (secondary, i.e., creep settlement) consolidation. Field

settlement and pore pressure data are compared with numerical predictions. Chapter 11

contains a summary of the work that has been presented, as well as conclusions and

recommendations for future investigations.











CHAPTER 2
LITERATURE REVIEW

2.1 Consolidation

The first author to deal with the important problem of fluid-filled deformable

porous solids was von Terzaghi. In a famous paper presented to the Academy of Sciences

in Vienna in June 1923, von Terzaghi showed the derivation of his consolidation theory.

This theory was later published in his book, which is now considered as the first

substantial book [30] in soil mechanics.

In the early 40's Biot [1] generalized von Terzaghi's theory of consolidation by

extending it to the three-dimensional case and by establishing equations valid for any

arbitrary load varied with time. In the following years, Biot generalized his theory to

include properties of anisotropy, variable permeability, linear viscoelasticity, and the

propagation of elastic waves in a fluid saturated porous solid [2,31]. The main

disadvantage of Biot's model, however, lies in the fact that the corresponding theory is not

developed from the fundamental axioms and principles of mechanics and thermodynamics.

Thus, some derivations are complicated and obscure. Finally, Biot [32] developed, within

the framework of quasi-static and isothermal deformations, a theory of large deformations

of porous media.

Since the beginning of the 1960's the study of porous media advanced with

development of new continuum theory of mixtures. In 1960, Truesdell and Toupin [26]

presented a treatise on the classical field theories where they developed in detail the








properties of motion and fundamental physical principle of balance. In 1965, Green and

Naghdi [25] developed a dynamic theory for the relative flow of the two continue based

on an energy equation and an entropy production inequality for the entire continuum.

With the advances in modem computational science and the development of

rigorous numerical techniques, such as the finite element method, numerical

implementations of the consolidation theory, Biot's equations and the mixture theories

found wide applications. A variational formulation of the dynamics of fluid-saturated

porous solids was the basis of a numerical method that Ghaboussi and Dikmen [33]

developed for the purpose of discretizing the spatial media into finite elements. Sandhu

and Wilson [34] first applied the finite element method to study fluid flow in saturated

porous media. With the presence of finite element method as a sound numerical technique,

it was possible to extend the mixture theory to encompass elasto-plastic non-linear

constitutive models and obtain reliable solutions of the field displacements and pressures.

A general analytical procedure that accounts for non-linear effects was presented by

Prevost [8]. In his work, Prevost focused on the integration of the discretized field

equations based on the mixture theories of Green and Naghdi [25]. Later he worked on

several numerical applications to study the consolidation of inelastic porous media [35]

and the non-linear transient phenomenon [36]. Due to the increasing necessity of nonlinear

applications, Zienkiewicz and other researchers published a series of papers that elucidated

various numerical solutions for pore-fluid interaction analysis. Zienldewicz et al. [37]

classified different method of analysis in a comprehensive paper on numerical solutions of

the Biot formulation. A continuum theory of saturated porous media that is applicable for








soils exhibiting large strains was formulated later by Kiousis and Voyadjis [38]. Borja and

AlarCon [39] recently proposed a framework for large strain consolidation based on

continuum theory of mixtures.

The mathematical basis for balance principles presented in this study is derived

from the general theory of mixtures [21-26]. This research is focused on fully saturated

porous media (two-phase continue). Field equations governing the interaction between

soil skeleton and pore fluid are developed from balance laws.



2.2 Large Strain Multiplicative Plasticity

Up to the beginning of the 1980s computational methods for large strain

elastoplasticity typically relied on hypoelastic extensions of the classical small strain

models; see, e.g., the reviews of Needleman and Tvergvaard [40], hence remained

restricted to small elastic strains. Computational approaches based on the multiplicative

decomposition appear to have been first proposed by Argyris and Doltsinis [41] within the

context of so-called natural formulation. Subsequently, however, these authors appear to

favor hypoelastic rate models on the basis that multiplicative formulations '.... Lead in

principle to non-symmetric relations between stress rates and strain rates' (see [21, p.

22]). Simo and Ortiz [43] and Simo [44] proposed a computational approach entirely

based on multiplicative decomposition and pointed out the role of intermediate

configuration in a definition of the trial state via mere function evaluation of hyper-elastic

stress-strain relations. Extensions of classical volume/displacement mixed methods within

the framework of the multiplicative decomposition, originally introduced for plasticity








problems in [45] are presented in [46]. In recent years, computational approaches based

on multiplicative decomposition have received considerable attention in the literature.

Simo [47] exploited a strain-space version of the principle of maximum dissipation to

obtain the associative flow rule consistent with multiplicative decomposition, and used a

(covariant) backward method to derive a large strain version of the return mapping

algorithms. Subsequently, Weber and Anand [48] and Eterovich and Bathe [49] used the

multiplicative decomposition in conjunction with the logarithmic stored energy function

and an exponential approximation to the flow rule cast in terms of full plastic deformation

gradient. The multiplicative decomposition along with a logarithmic stored energy

function is used in [50]. More recently, Moran et al. [51] addressed a number of

computational aspects of multiplicative plasticity and presented explicit/implicit integration

algorithms. Methods of convex analysis, again in the context of the multiplicative

decomposition, are discussed in [52]

The preceding survey, although by no means comprehensive, conveys the

popularity gained in recent years by computational elasto-plasticity based on the

multiplicative decomposition. Despite their success, these approaches involve

modifications, and often a complete reformulation, of the standard closest-point

algorithms of the small strain theory. From a practical stand-point the implication is that

the implementation of classical models needs to be considered on a case-by-case basis in

the large strain regime.

In a later study, Simo [20] proposed a state-of-the-art algorithm based on

multiplicative decomposition of the deformation gradient, as suggested by Lee [53],








Mandel [15] and others. The appeal of his formulation is that: the closest-point projection

algorithm of any classical simple-surface or multi-surface model of small strain plasticity

carries over to the large deformation context without modification. In particular, the

algorithmic tangent moduli of small strain theory remain unchanged while introducing a

further simplification: the closest-point projection algorithm is now formulated in principal

(Kirchhoff) stresses. For the static problem, the proposed algorithm preserves exactly

plastic volume changes if the yield criterion is pressure insensitive. For dynamic problems,

Simo [20] presented a class of time-stepping algorithms which inherits exactly the

conservation laws of total linear and angular momentum. Simo's method is adopted with

success in some recent works [39,54,55] for formulating nonlinear plasticity model in

large deformation context.

Present study has followed the above-mentioned theory [20] of multiplicative

decomposition to derive explicit expressions of consistent tangent moduli for the proposed

elasto-plastic and viscoplastic constitutive models (see Chapters 5 and 6).










CHAPTER 3
THEORY OF MIXTURES


3.1 Basic Theory

The conceptual model of a multiphase continue is based on the phenomenological

behavior of each phase rather than particulate nature and the microscopic origin of the

phenomenon involved. In other words, each phase (or constituent) enters through its

average properties obtained as if the particles were smeared out in space. In order to be

able to derive multiphase field and constitutive equations for such a medium, a technique

for obtaining local average quantities is therefore necessary. Furthermore the basic

kinematics and balanced equations for each phase and for the mixture as a whole must be

defined. In following sections (3.2 to 3.4) general kinematics and balance laws are derived

for a multiphase (n-phase) continuum allowing for the selection of the constitutive

relations to be defined according to the particular phases that composes the mixture. For

simplicity, it is assumed that phases are non-interacting and non-polar.



3.2 Kinematics

A mixture can be viewed as a superposition ofn single materials each of which

may be regarded as a continuum. It is assumed that at any time t each place x of a mixture

is occupied simultaneously by n different particles: X', X2, ...., X". As in single-phase








theory, a fixed but otherwise arbitrary reference configuration and a motion are assigned

to each phase [26] as


S= (Xa)


V a=1 ...n,


Figure 3.1 Geometric representation of kinematics for a two-phase mixture



where X" denotes the position of a phase in its reference configuration, and x is the spatial

position occupied at the time t by the particle labeled X. The function 4a in (3.1) is


called the deformation function for a phase at time t. In classical continuum theory 4a is

assumed invertible, and thus


(3.1)









X = a (x) V a=1...n. (3.2)

The invertibility of the deformation function ensures that a particle at X" cannot occupy

two spatial positions at a given time and that two particles of a phase with positions

X and X cannot occupy the same spatial position. Figure 3.1 shows geometrical

representation of (3.1) in Cartesian coordinate system e R3.

The velocity and acceleration of X" at time t are obtained from Equation (3.1) by

time differentiation, viz.
Da
vi = v(X, t) ( i)i;
t (3.3)
a = a(x, t)= D i) =
Dt

where a superimposed dot indicates differentiation with respect to time holding X" fixed

(i.e., the material derivative following the motion of a phase). Material time derivative

may be computed from spatial description using the following definition

Da )= )+ va gradee. (3.4a)
Dt at

Material time derivative of a volume integral can be expressed as

DG () adn= (*) dl + ()a (a n) dr. (3.4b)


See [56] for proof of the identities (3.4a), (3.4b). Here, and in the following, grad and

GRAD are used to denote spatial and material derivatives, respectively. The deformation

gradient for X" at time t is defined by









F = Fa(x,t)= GRAD = a (3.5)
.A-

and the velocity gradient is defined by

a = la(x, t)= gradva =Fa (Fa) (3.6)



3.3 Average Quantities

An important assumption in the theory of mixtures is that the phases of a mixture

are allowed to occupy common portions of a physical space. Then each spatial position x

in a mixture is occupied by n elements, one from each phase (see Figure 3.1 for n = 2). To

address this assumption one needs to define average quantities.

Average quantities are obtained by integrating microscopic quantities over an

averaging volume or area. In the macroscopic field, the averaging volume represents a

physical point, denoted by dV. Similarly, the averaging area dA, represents and

characterizes a physical point on the surface of dV, and is an infinitesimal element of area

in the macroscopic field. The part of dV occupied by the a phase is denoted by dV", and

the volume fraction n" of the a phase is the fraction of dV occupied by the a phase

defined by


na =n(x,t)= dV (3.7)
dV

where n" is constrained by X n' = 1 and 0O n" 5 1. Similarly, the part ofdA lying in the
a


a phase is denoted by dAY, and the real fraction is defined by









a a dAa
n = n (x,t) (3.8)
dA

Again, Z W" = 1 and 0 _< a <_ 1. It is assumed in the following that the identity
a

o" = no holds [57].

A macroscopic average mass density function, p" is associated with each phase and

is defined as the volume average of the microscopic density function, pa.

P dV1 Pdu, (3.9)


where du is the microscopic volume element. The intrinsic volume average mass density is

defined as

-ia 1 I adu p". (3.10)
dVadV n

If the mass density of the a phase is microscopically constant the intrinsic volume average

mass density equals to the microscopic mass density, i.e., pa = pa and thus pa = nepa.

The mass density p of the mixture is defined as

p=p(xt)= p, (3.11)




and the volume-average velocity v for the mixture is defined as

S= i(x,t)= pava. (3.12)
Pa








Following the similar averaging approach, a macroscopic partial stress vector t" may be

defined as

ta =1 L tada, (3.13)
dA dA

where da is microscopic area element. t. denotes the intrinsic stress vector of a phase,

t" = n"t .



3.4 Balance Laws

Once kinematics and local average quantities of a mixture are derived, one may

postulate the laws of balance based on the theory of mixture [26] which must be satisfied

irrespective of constitutive relations. Phases are understood to be material elements, which

are open systems on a local state. Accordingly, local balance relations are derived for each

individual phase. The equations are obtained in spatial configuration by applying the

fundamental laws of mechanics: balance of mass, balance of linear momentum, balance of

angular momentum, and the first and second laws of thermodynamics.

For consistency in notations, in the derivation of balance laws and in the following

thereafter, spatial (deformed) configuration is represented by domain R, bounded by

surface I while material (undeformed) configuration is represented by domain B, bounded

by surface B.








3.4.1 Balance of Mass

The balance of mass can be expressed as "The time rate of mass in a fixed region in

space 0 is equal to the time rate of mass flowing through the surface r that encloses G."

In equation form

J dn = -Jrpv ndr, (3.14)

where n is the unit outward normal to the surface F, v is the flow velocity.

For illustration, mass flow through a control cube (Sxx8yx5z) in Cartesian

coordinate space e R' is considered in Figure 3.2. Surface areas normal to the flow

components along coordinate axes are Fr= 6ysz, Fy = 6x8z, z = 85xy.




Pa Vy Ay" no PaVzn-zna




PaVx-Axn PnaV.V+Ax.nO
X

z
PaVz+Azn ----


P. y-,Ayn


Figure 3.2 Balance of mass: flow through a control volume








For the cube (3.15) can be written as,

S a a Pavo a sySz + av+Ay na -PaVy a6x
Ax+Ax PaVx-Axn y+Ay -paVyA(3.15)
+ av na _pXva_ zna]8x~y= pdV t+,t P adVu t
+Az -p -Az =t+At it

Dividing by the volume dV = 8xSy8z and At, then upon rearrangement (3.16) can be

written as

div(naPav a+ anaPa) = 0. (3.16)

(3.16) is the localized form of balance of mass for a phase.



3.4.2 Balance of Linear Momentum

In order to establish the balance of momentum laws for each a phase, one needs to

consider the forces acting on the a phase within the region 0 such as drag forces, body

forces or gravity forces, as well as, the effect on the a phase of the mixture outside the

region f. This effect is accounted for by introducing partial stress vector t, defined in

(3.13). Now, let 0a be the Cauchy partial stress tensors [21-26] and n denotes the unit

normal vector to the surface F. ao is related to t* by the relation t' = a n.

The first Euler equation postulates that "The rate of change of the total momentum

of a given mass is equal to the vector sum of all the external forces acting on the mass." In

equation form


U Ipavadf= Jfadn+ Jf a -ndr+ fpagdo. (3.17)
Dt fa*n*r o









g is the vector of gravity acceleration, p" is the momentum supply for the a phase from

the rest of the mixture due to other interaction effects, e.g., relative motion of the phases.

pa is subject to a = 0.


From divergence theorem

f a n dI'r= div aodn. (3.18)

Substitution of (3.18) in (3.17) results in


D p Lvad = pa D- va)dn=J[divaa C ++p Pag d, (3.19)

which simplifies to the local version of the balance of linear momentum equations for

individual phase as

paaa = divoa +a +pag. (3.20)



3.4.3 Balance of Angular Momentum

The principle of balance of moment of momentum states that the material rate of

change of moment of momentum of a body about a fixed point xo, is equal to the resultant

moment acting on the body around that point.

The moment of momentum for a phase is defined as

La = J Opr x vadl. (3.21)

r = vector of moment arm = x xo. Making use of (3.17) material derivative of (3.21) can

be written in the form








DaL Da Ip r x vad
Dt Dt (3.22)
= rx (pg+a dn +frrx (o .n)dr.

From divergence theorem

I rx ( a .n)dF = Ir nrxa c)dF =I div(rx aa)dn. (3.23)

Substitution of (3.23) in (3.22) yields

D L Pa rPa rx vda
Dt Dt J(3.24)
= J{r x (Pag+ +div(r xaa)df.

Local form of (3.24) results

pa Darxv' =rx (pag+ a) + div(rx c1. (3.25)

For cases of non-polar phases partial stress tenor, a" are symmetric since there is no

supply of moment of momentum (i.e., antisymmetric part of o is zero).


3.4.4 Balance of Energy
Principle of balance of energy or the first law of thermodynamics states that the
material rate of change of internal energy in a body is equal to the resultant deformation
power acting in the body plus the rate of heat added to the body.
The internal energy in the body for a phase may be defined as

Ua = Jpae 0 d, (3.26)








where e" represents internal energy density for a phase. Assuming no supply of energy

due to interaction, e.g., chemical interaction within the phases material rate of change of

If can be decomposed into following components

Du Dt P +K( +Qao (3.27)
Dt Dt

where P" = Mechanical energy due to deformation of a phase

K" = Kinetic energy of a phase

Q" = Thermal energy of a phase due to heat generation within the domain 0 and

heat flow across the boundary dF.

Material rates of energies for a phase can be expressed as follows:

Dt P pag. vadn+Ja .vadn+ frOa va @ndr, (3.28a)


DaK Da 1 a a
Dt Dtj 2 (3.28b)
a-fn lP va.a d + r (lpava.v va .n)dr,



Dat CaHa d-f r qa.ndr (3.28c)


If is the heat generation per unit mass for the a phase, and qC represents the heat flux

vector associated with each phase. Substituting (3.28) in (3.27), one obtains for a phase

the following expression








De Ivpa eadn= fpag.vcd+J I a.adn+Jro :vaC ndr

+-a I P a v a d.v + r ava .v~v v a -n)dr (3.29)

+ IJ aH d q-Irqa -ndr.



3.4.5 Entropy Inequality

Entropy inequality or the second law of thermodynamics puts limit on the direction of

such processes where thermal phenomenon are involved, permitting energy transfer to occur

spontaneously only in certain preferred directions. The limitation may be expressed

mathematically as an inequality stating that the intrinsic entropy production of the entire

mixture is always nonnegative and is positive for an irreversible process. In other words the

material rate of change of the entropy increase is always higher or to that of the entropy due

to heat transfer.

The entropy density for the entire mixture, it can be defined as

11 = EP 1a (3.30)
a

TI" is the entropy density for a phase. Entropy increase for the entire mixture results

JpTidd= fI pa ad dn. (3.31)
a

Assigning to each phase a temperature 0, given by a positive-valued function, the second

law of thermodynamics may be expressed as


p t d- E d + J 1 -q ndr 0, (3.32)
a a a a O






24

FH, q" are defined in (3.28c). Using (3.4b), divergence theorem and balance of mass

(3.16), material time derivative of jrparl df can be derived as


L paLdn = pa a+div[Pa ava d


= {ak Pa +div pav) +Pa( Tla +gradna .va} dn (3.33)

at


Substitution of (3.33) in (3.32) results


a a p + div1 q dQ 0 (3.34)
a t Oa" .00C

Localization of (3.34) yields

pa Nt ae + div a ]>0. (3.35)
at a Iw[ 8a J


Introducing the Helmholtz free energy function for a phase, va defined by

acc ea oaCr1, (3.36)

one can write

p a a 1&.a / & 01 SB 01 y
p a Z t (3.37)
at Oa 9t tt

Substituting (3.37) in (3.35), localized form of entropy inequality can be written as

pa a ~a a JPaH a +adiv 1qat 20. (3.38)
at a at O








3.5 Balance Laws for Saturated Soil

Balance laws of saturated soil medium are special case of generalized balance laws

of n-phase continuum considering n = 2 (solid and fluid phases). Following the derivations

in Section 3.4, necessary balance laws for soil-water mixture are deduced in this section.

Motions of solid and fluid phases are considered separately.



3.5.1 Balance of Mass

Following (3.16), localized form of balance of mass for solid and fluid phases can

be written, respectively, as

div(nspsvs)+ (nsps) = 0, (3.39a)


div(nWpwvw + (nwpw = 0. (3.39b)

Adding (3.39a) and (3.39b) gives the conservation of mass for the soil-water mixture as

p+div(p)= 0, (3.40)
at

where i is the volume-average velocity, defined in (3.12).



3.5.2 Balance of Linear Momentum

In the absence of inertia forces balance laws of linear momentum for solid and fluid

phases can be written from (3.17) as follows:

Jpsg df+J s d++J fc-ndrf= (3.41a)


Jpgd +r w d+r owT nd=0. (3.41b)









Since Os and Ow are the internal forces which naturally will not affect the soil-

water mixture as a whole, ps + w = 0, i.e., the seepage force exerted by the fluid on the

solid matrix is the negative of the reactive force exerted by the solid matrix on the fluid.

Consequently, the sum of (3.41a) and (3.41b) results in

Spg dl + nd = 0. (3.42)

o is the Cauchy total stress tensor, a = s + .

Now, let P and P' be the (non-symmetric) first Piola-Kirchhoff partial stress

tensors arising from the fluid and intergranular stresses, respectively. Also, let N denote

the unit normal vector to the surface aB of the undeformed region B. The tensor P" is

defined such that P" N represents the resultant force exerted by the fluid per unit area of

the solid matrix in the undeformed configuration. Similarly, the product P' N is the

resultant net force exerted by the individual grains (which may include the partial effects of

fluid pressures) over the same undeformed reference area. By the additive decomposition

of the Cauchy total stress tensors, we obtain a similar expression for the first Piola-

Kirchhofftotal stress tensor P

P= J. F-t = p +Pw, (3.43)

where Ps = Joa F-t and P" = Joa F-t are the first Piola-Kirchhoffpartial stress

tensors arising from the solid and fluid stresses, respectively, and

J = det(F); F= x=X+u. (3.44)
dX







In (3.44), J is the Jacobian, F is deformation gradient, X are the coordinates of a

point X in undeformed configuration, u is the macroscopic displacement field of the solid

phase, x are the spatial coordinates of point X.

P can also be decomposed as

=P +- (3.45)


where P is the first Piola-Kirchhoff effective stress tensor, and (P" /n")-N represents the

resultant force exerted by the fluid per unit area of void in the undeformed configuration.

P and P" are related by the equation

Ps =P+ 1 Pw. (3.46)


An integral equation similar to (3.42) can be developed in terms of the tensor P.

With reference to the undeformed configuration, (3.42) can be written in the form

JBPogdV + fB .NdA= 0, (3.47)

po = Jp is a pull-back mass density of the soil-water mixture and g is the vector of gravity

accelerations.

It is very important to note that po is not a constant quantity. Figure 3.3 would

explain the scenario. In Figure 3.3, the fluid now occupying the void in a soil at a point

4(X) may not necessarily be the same fluid material that occupied the same void at a

reference point X in the undeformed region B. Mathematically, t(X) = t(Y), where 4t(Y)

is the spatial configuration of the fluid phase whose reference configuration is at Y.











dV, ddY





YY)




JdV, dr'





Figure 3.3 Balance of mass: total mass of material point X is not conserved in (t(X)



Notice that Y does not necessarily have to be in B. Likewise, fluid phase of material point

X, X' might not necessarily be present in the spatial configuration 0t(X). Hence, po does

not necessarily represents the true mass-density in B of the soil mass which now occupies

the volume 4k(B), since fluid can migrate into or out of the soil matrix during the motion.

In other words, the total mass of the soil-water mixture in B is not necessarily conserved

in C(B).

A simple relationship analysis in the following would demonstrate the effects of

diffusion on mass densities. Let n1 (X, t = 0) be the initial porosity of the point X in B.


Then, the initial volume of the voids in an elementary volume dV is no dV, while the









initial volume of the solid region is ( n' )dV. As the soil matrix deforms, its volume

changes to dv = J dV. Now, assume the solid phase is incompressible. Since u is the

displacement field of the solid phase, its volume is conserved at (1- nw dV in dv, while

the volume of the voids changes to dv (1 n )dV. Consequently the porosity varies

according to


nw 0 -d- (-(nw) -1. (3.48)
JdV 0

Hence, the total mass density and the porosity of soil vary with deformation through the

Jacobian J.


3.5.3 Balance of Energy

Ignoring kinetic energy and non-mechanical power, and assuming balance of

momentum and balance of mass hold, (3.29) can be simplified for balance of energy of

a phase as


fo paedn= Jpag.v adO+ Ja vadn+ f ra va ndr. (3.49)

The localized version of balance of energy can be derived in the following fashion.

Consider the left-hand side of (3.49), for example. Following the derivation of (3.33), one

can reduce

D I paeadn = jpaea dn, (3.50)








since grad e = 0. Using divergence theorem and localized version of balance of linear

momentum (3.20), one can reduce right-hand side of (3.49) as


J(pag+a-a)vadn+ fro : a ndr
-odivao vadn+fo [div +ava+o :gradv ]d (3.51)

= j a : a dQ= Jca :da da

I" is defined in (3.6), d" = Sym(l"). Since a" is symmetric, a" : I' = a : d". Substituting

(3.50) and (3.51) in (3.49), localized version of balance of energy for a phase can.be

written as

paa = oa :da, (3.52)

Corresponding localization for the saturated soil media takes the form

p = :d' : + W :dW, (3.53)

where e is the rate of internal energy for the soil-water mixture obtained from the volume

average

pSOS pWiW
e= s (3.54)
P

It is often convenient to describe the balance of energy in the material picture

because the domain of integration of the functions remains fixed. To this end, one makes

use of the following transformation. Let the right leg of the tensorP be pushed forward by

the configuration 4i. The result is the Kirchhoff total stress tensor -, which differs from

the Cauchy total stress tensor by the factor J, i.e.,

r=J= P]Ft. (3.55)









T can be decomposed into solid and fluid counterparts in any of the following ways

,w
T= s + tw = +-. = +1. (3.56)
nw

r is Kirchhoff effective stress, 0 is Kirchhoffpore water pressure (compression positive), 1

is the second order identity tensor. Now, pulling back the left leg of P by the inverse

motion (,t)'1, one obtains symmetric S, the second Piola-Kirchhoff total stress tensor

such as

S= F- -P = F- -F-t = JF-1' -iF-t. (3.57)

Following (3.56) and (3.57), additive decomposition of S can be written as

S= S+9C-1, (3.58)

where S is the second Piola-Kirchhoff effective stress tensor and C is the right Cauchy-

Green tensor given explicitly by

C= Ft F. (3.59)

F is defined in (3.44).

Let x = |(X, t), E'(X, t) = eL(x, t). If one multiplies the localized balance of energy

(3.52) by J, and uses the porosity expression (3.48), one obtains the following expression

for the balance of energy for the solid and fluid phases in localized material form

PsnsES = :ds = Ss :t;
(3.60)
pw J-n )tw =rw= :dw SW .
0 2








Balance of energy for the soil-water mixture in the material picture is now given by

JpE= r8 :d' +9w :d" = S' :C+-S"W :CW, (3.61)
2 2

Where E is obtained from the volume average


E +Ps( 0 e. (3.62)
Jp

The quantity Jp E is the mechanical power generated per unit reference volume of the

soil-water mixture. (3.61) can be expressed in a more elegant form in terms of effective

Kirchhoff stress and deformation of the solid phase as

JpE=': d', (3.63)

i.e., the sum of the mechanical powers of the partial stresses is equal to the mechanical

power of the effective stresses with respect to the deformation of the solid matrix

computed from its own motion.

Proof of (3.63) is given in section A.2. (3.63) states that total mechanical power in

soil-water mixture is absorbed by the energy rate : d', and that the pore pressure tensor

x" / n" in (3.56) performs no work. It is obvious from the fact that fluid is assumed

incompressible and has no shear strength. By virtue of these assumptions, fluid cannot

store volumetric nor deviatoric energy, i.e., it has no mechanical power.



3.5.4 Entropy Inequality

Ignoring non-mechanical power and kinetic energy production, the localized form

of entropy inequality or the second law of thermodynamics (3.38) can be simplified as









pa fa e0 l 20. (3.64)

In case of soil-water mixture (3.64) takes the form

e- 7 0. (3.65)

e is defined in (3.54). Similarly, J is the rate of free energy for the soil-water mixture

defined as

pS S W+p W
S= (3.66)
P

Let x = >(X, t), Y"(X, t) = w"(x, t). Similar to (3.62), volume average rate of free energy

T can be expressed as

psns s +ps J-nw )lw
JY = (3.67)
Jp

From (3.62), (3.65), and (3.67) localized version of entropy inequality for soil-water

mixture in undeformed configuration B can be written as

JpE JpY 2 0. (3.68)

JpY is the power generated from free energy. JpY = d Y /dt; Y denotes the stored

energy function, or free energy, per unit reference volume of soil matrix. Substituting

(3.63), one can rewrite (3.68) as

xr:ds 20. (3.69)
dt









CHAPTER 4
VARIATIONAL EQUATIONS, CONSTITUTIVE THEORIES



4.1 Boundary Value Problem

In order to formulate a well-defined boundary value problem for consolidation

phenomenon, one needs to consider a problem domain with a set of suitable boundary

conditions. Let B c Rmd (nsd = no. of spatial dimensions) be an open set in material

configuration (time to) with piecewise smooth boundary 8B. MB is assumed to admit the

following decomposition

aBdUaBt=aB, t.
{ % B =, (4.la)
aBd n aB = 0,
SuaBb = aB (4. b)
s1e n Bh = 0.

Bd', t8, sB9, 0Bh are open sets in 8B. sBd and 9BB represent the portions of 9B with

prescribed displacement and traction, respectively while oB6 and oB" represent portions

with prescribed fluid potential and volumetric flow rate, respectively. 0 is a null set.

In spatial description at any time te[tn, t.n+], reference domain B will be mapped

to the configuration 0 = t(B) c R"d with boundary F = t(8B). Similar to (4.1),

decomposition of boundary F will now take the following form

Srd rt = (4.2a)
d nrt =0,









riT Fh = 0
r0 nph = 0.


(4.2b)


-


Figure 4.1 Prescribed boundary conditions of spatial domain C


In (4.2a) and (4.2b), I" = t(aBd), rI = (aSB), Fe = 4(aBe), rh = t(aBh). Figure 4.1

shows decomposition of domain boundary in spatial configuration. Having outlined the

domain boundary, one would require strong form of consolidation phenomenon with

appropriate boundary conditions to define a well-posed mathematical problem.


4.2 Strong Form

Strong form or field equations of consolidation problem emanate from balance

laws. Equation of equilibrium is derived from balance of linear momentum. Localization

of (3.47) results in the following statement of stress equilibrium in material description








DIVP + pog = 0. (4.3)

DIV is the material divergence operator. (4.3) is subject to the following boundary

conditions: the motion 0 is prescribed to be Od on aBd C aB, and the traction P N = t is

prescribed on the remainder aBt; and subject further to the constraint imposed by the

balance of mass. Push-forward of (4.3) in spatial configuration E can be obtained as

div i + Jpg = 0, (4.4)

since grad J = 0 (see Section A. 1).

Equation of flow continuity is derived from the balance of mass. It is assumed in

this study that both the fluid and solid phases are homogeneous and incompressible, i.e.,

a(pa)/at = 0; grad pa = 0 (see Section A.1). Using these assumptions, pa can be factored

out and eliminated from (3.16), resulting an expression

divnC )+ (nO)= 0. (4.5)

Adding (4.5) over the phases yields

div i- nW)v +div(WvW)= 0. (4.6)

For future reference, it is useful to define a superficial, or Darcy, velocity as

= nw(vw -s). (4.7)

The vector V represents the relative volumetric rate of flow of fluid per unit area of

deforming soil mass in spatial configuration 0 = 4((B). V is related to fluid potential H

by constitutive relationship. See Section 4.5 for discussion. For simplicity of notation, v

will be used for solid phase velocity, v' in subsequent discussion. Substituting (4.7) in

(4.6), field equation of flow continuity can be obtained in the spatial reference as








div v + div v = 0. (4.8)

(4.8) is subject to the following boundary conditions: fluid potential H[ is prescribed to be

ni c: t(aBB), and the volumetric flow is V n = -q on the remainder *t(BBh); and subject

further to the constraint imposed by the conservation of momentum. Here n is the

outward unit normal to the deformed surface 0(#B) and q is positive when fluid is being

supplied to the system. In material description, (4.8) takes the following from

DIV V + DIVV =0. (4.9)

V, V are pull-back velocities in undeformed, reference configuration B. V, V can be

obtained through Piola transform of v, V such that, V = JF~- V, V = JF-' v. Let

V N = -Q be the prescribed volumetric rate of flow per unit undeformed area across the

boundary OBh, N being outward unit normal to the undeformed surface aB. Q maintains

the same sense of direction as q.



4.3 Weak Form

In order to establish weak form of the boundary value problem, one needs to

define following spaces in accordance with the standard arguments of variational

principles. Let the space of admissible configurations be

C4 = {O: B Rnsd 4i eH'1, = Od on Bd (4.10)

and the space of admissible variations be

V4 = r:t B -> R nsd i eH1, T= 0 onaBd}, (4.11)








where 'r is an admissible virtual displacement field, H1 is the usual Sobolev space of

functions of degree one. Further, let G : C4 x V4 R denote the weak form of

equilibrium equation (4.8) in material description.

G(,H,nT)= fB(GRADTI: P-poTig)dV- Jg.t dA. (4.12)

The balance of linear momentum is given by the condition G($, I, 1) = 0. The formal

statement of (4.12) is: Find 0 e C# such that G(4, H, 1) = 0 for all ie V4. Using (3.55)

and (3.56), one can rewrite an equivalent expression of G of(4.12) in the same

undeformed, material configuration with the integrands evaluated in spatial description as

follows:

G(0, H, ) = Js(grad T :r + div i-poT. g)dV- Jls T t dA. (4.13)

Now, let the space of potentials in spatial reference be

Ce = H: t(B)-+ Rs neH1, n= ig onrO} (4.14)

and the corresponding space of variations be

V = Vw:t(B)>Rnsd eH'I, w= onrO}, (4.15)

where W represents an arbitrary virtual pore pressure field. Further, let H : Ce x Ve-+ R

denotes the weak form of equation of continuity (4.8) in spatial description.

H(,H,w) = f(\wdivv -gradyr- )df- lJw- qdr. (4.16)

One can show that the balance of mass is given by the condition H(, H, v) = 0. Formal

statement of weak form (4.16) will translate as: Find II Ce such that H(, I, ) = 0 for

all admissible We Ve.








The domain of integration of H(<, H, i) of (4.16) can be of evaluated quite easily

in undeformed, material configuration by introducing the Jacobian J. In doing so, one

would require the following identity

J= Jdivv. (4.17)

J is the time derivative of the J. See Section A.4 for proof of the identity (4.17).

Substituting (4.17) in (4.16) yields

H(,ln,w)= J(vj gradvy -) dV- J wQ dA. (4.18)

Relation of area transformation of flow rate, q dr = Q dA (see Section A.4) is employed

in deriving (4.18). H(<, H, W) in material description takes the form

H($,H n,)= JB (- GRAD V) dV b vQ dA. (4.19)

(4.18) and (4.19) are equivalent expressions since GRAD V = grad N.

Presence of rate term J in the variational equation H(, n, ) makes it

mathematically awkward. One can eliminate rate terms altogether by semi-discretization

in time domain, via finite difference, for example. Following is a brief discussion on time

descretization scheme for consolidation problems.



4.4 Time Descretization Scheme

The ordinary differential equation associated with the problem of consolidation is

generally stiff. A physical insight provides an explanation: points near drainage

boundaries consolidate many times faster than do points at remote places. The spectrum

of eigenvalues associated with the consolidation equation is therefore wide.








The general linear k-step method for approximating the solution of a system of

ordinary differential equations of the first order,

a=f(d,t), d(0)=do, t>0 (4.20)

is

k
(acmdn+i-m + Ati3mfn+l-m)= 0, (4.21)
m=0

where At = t,+ -t, and am, 3. are unknown coefficients. A linear multistep (LMS) method

is explicit if Po = 0, otherwise, it is implicit. Implicit method is preferred in this study

because it has a larger region of stability than the explicit methods, and it is compatible

with the stress point algorithms used in the development of constitutive models (see

Chapters 5 and 6). As a result, 1o 0. (4.21) has 2k + 2 unknown coefficients.

An effective technique for solving stiff differential equation is provided by so-

called stiffly stable methods proposed by Gear [58-61]. These k-step methods of order k

are based on backward differentiation formulas (BDF) derived by setting bo 0, 01 = 32

=... Pk = 0. The resulting BDF approximation is

k
(amdn+l-m)+ AtPofn+l = 0 (4.22)
m=0

with k + 1 unknown coefficients which one can choose by forcing a k-step method to

satisfy an accuracy of order k. There is an arbitrary normalizing factor so that one can set

to = -1.0, leaving k unknowns. (4.22) then takes the form


dn+l (amdn+l-m) -Atlofn+l =0. (4.23)
m=1 ,









See [5] for determination of coefficients for BDF-k scheme. (4.23) can be made more

generalized by introducing one-step recurrence relation for f such that


dn+i- ((amdn+i-m) -AtPo[lfn+i+(1-P)fn = 0. (4.24)


O < p 1 is the parameter of generalized trapezoidal method. Some of the well-known

families of 1-methods are presented in the following.



Table 4.1 Generalized trapezoidal method

P Method
0 Forward Euler
1/2 Crank-Nicolson
1 Backward Difference



Unconditional stability is achieved for any At if 1 > 1/2. In general, the nonlinear

responses of interest are dominated by low-frequency component of the system, but high

frequencies also enter into the solution because of the numerical approximation. It is

known that the backward difference scheme (P = 1.0) can damp such high-frequency

components but at the expense of accuracy. On the other hand, Crank-Nicholson scheme

(1 = 1/2) possesses a second order accuracy but lacks the numerical dissipation of the

backward difference scheme. It was shown in [61] that the variable step size, variable

order BDF methods are convergent and unconditionally stable for ordinary differential

equations.

Following (4.24), one can obtain time-integrated variational forms of

HAt (, I, w) given as









HAt Hn.W) =v/ -f Jn+i- SamJn+i-m dV

-P0 oJb(grad JV)n+ +(1 )(gradl .J-)n] dV (4.25)
-Po10BPQn+l +(l-P)Qn]dA= 0,


HAt(,Hw)= Jn+l mamJn+-m dV

-PoJg0 p(GRADyV)++(1-P)(GRADw.-)n]dV (4.26)
-Poj oIPQn+ +(1-P)Qn]dA=0.

(4.25) and (4.26) are obtained from (4.18) and (4.19), respectively.


4.5 Large Deformation Plasticity Model for Soil Skeleton

In this study, plasticity behavior of soil skeleton in large deformation is based on

multiplicative decomposition of the deformation gradient, F [15, 16]. Let X be a

macroscopic point containing a sufficient number of solid particles in the reference,

undeformed configuration B, and x be the configuration ofX at some time t > 0, i.e., x =

t(X). Recall from (3.44), F = Ox/8X. x and X are coordinates ofx and X, respectively.

The motion of X produces both reversible as well as irreversible microstructural

changes in the soil. Typical processes associated with reversible microstructural changes

include elastic deformation and (for plate-like particles) elastic bending of the granules

comprising the assembly. As x is unloaded, it moves to some intermediate, stress free

configuration defined by the macroscopic point x". Assuming that this intermediate

configuration exists, the chain rule can be used to express F in the product form










F=- u Fe. FP


V = ox/Ox", P = 8x"/aX. Figure 4.2 presents the schematic of the multiplicative

decomposition ofF.


Reference
Configuration

X


(4.27)


x = (X)
Spatial
Configuration


Intermediate
Configuration


xe

Figure 4.2 Illustration of multiplicative decomposition of the deformation gradient



Ignoring non-mechanical power and kinetic energy production, balance of energy

and the use of the second law of thermodynamics lead to the following reduced

dissipation inequality

dY 1 d'Y
D = d- S C-dr 0 (4.28)
dt 2 dt









where d = Sym(l) is the rate of deformation tensor, I = F F-I is the overall spatial

velocity gradient, and W is the stored energy function. Clearly, dW/dt = JpE and D = 0

for an elastic material (see (3.63)).

The form of the stored energy function \ determines the constitutive

characteristics of the soil. For isothermal, elastic processes, y depends only on X and C if

it is to satisfy the axiom of material and frame indifference. Equally well one can say that

for isothermal, elastic processes W is a function of X and elastic left Cauchy-Green

tensor, b', provided that b* satisfies an objective transformation.

An elasto-plastic process requires a yield function, a hardening rule, and the

imposition of consistency condition. Let f be the yield function defined as f(z, ) = 0. X

e Rm is a suitable vector of m > 1 (stress-like) internal variables characterizing the

hardening response of the soil. t e R" is a vector of internal plastic variable (strain-like)

conjugate to X in the sense that X = 8 V/8 From the framework set in (4.27), takes

the form

y= b, );be; be= Fe.(e)t. (4.29)

Now, consider the following time-derivative of Y

dY d. e : db .: be+be. t+ be + S., (4.30)
dt dbe dbe 8(

where 1, be is the Lie derivative of b. Inserting (4.30) in (4.28) and assuming isotropy in

the sense that b and 8Y/l8b commute, the dissipation inequality can be expressed as

( dY d 1 1 .
D= -2-- bc :d +2- b: be -be- -2. (4.31)
dbe dbe ( 2 a2







The first term of (4.31) yields the constitutive relationship

dY
=2 d .bc, (4.32)
db6

while the other terms yield, following the requirement f(r, x) = 0 and the postulate of

maximum dissipation,
I e f e .f
-l be fe = -.b, = -. (4.33)
2 && ax

0 and f satisfy the requirements of consistency conditions such that p > 0, f < 0, and @ f

= 0. Thus (4.33) defines the flow rule. (4.32) and (4.33) satisfy the reduced dissipation

inequality of (4.28) even with the use of empirically derived hardening law. The flow rule

(4.33) possesses a number of important properties. In particular, it gives the correct

evolution of plastic volume changes as the following observations reveal:

(i) The total and elastic volume changes are given by J = det(F) > 0 and J' =

(det(b'))' > 0, respectively.

(ii) Let P = det(FP). The rate of plastic volume change predicted by the (4.33)i is

given by the evolution equation

d (log JP)= tr[ ], (4.34)

which implies exact conservation of plastic volume for pressure insensitive yield

conditions; i.e., iftr[Of/t] = 0.

The left Cauchy-Green tensor b0 can be decomposed spectrally into


be= A()2m(A), m(A) = n(A) n(A), (4.35)
A=1









where A is the elastic principle stretch corresponding to the principal direction n(A).

is a vector operator defined as (ab)ij = aibj for any vectors a and b. Since r and b"

commute, T can be decomposed spectrally in the form

3
S= .A(A), (4.36)
A=l

where PA are the principal Kirchhoff effective stresses.

By the assumption of frame indifference and isotropy, the free energy function

can be expressed as symmetric function of the elastic principal stretches, i.e.,

(Xbe= 1yX 2e3,, 3); i=ln(A), A=1,2,3, (4.37)

where se 's are principal elastic logarithmic stretches. Thus, the elastic constitutive

equation (4.32) reduces to a scalar relationship between PA and es such that


PA = ; A = 12,3, (4.38)


In the elasto-plastic regime, the additional task of enforcing the consistency condition,

f(t, X) = 0, is done incrementally. In the first step, plastic flow is frozen and an elastic

assumption ignoring the constraints imposed by yield criterion leads to elastic a trial

elastic state.

f=l.f; be = 2Sym(Ibe); -=0, (4.39)

where f = Ox/ix, is the deformation gradient evaluated relative to the converged

configuration n (B). In the second step, trial state is held fixed and plastic relaxation is

introduced. The algorithm is given explicitly









f=0; 6e = -2 be; = (4.40)

subject to @ 0, f 5 0, and 0 f = 0.

Incremental counterparts of the evolution equations (4.39) and (4.40) are obtained

from the so called product formula algorithm [62]. From (4.39), trial elastic left Cauchy-

Green tensor is obtained in incremental form by freezing plastic flow as

be,tr =fben t; =n, (4.41)


where be and tn are the respective values of b and t at configuration tn (B). Similar

to (4.27), b',r can be decomposed spectrally in the form


etr = f tr(A); tr(A) = tr(A) tr(A). (4.42)
A=1

Introducing the product formula algorithm into the plastic flow equation then yields

be =exp -2AP -b betr; = n +A(p (4.43)

where Ap is an incremental consistency parameter that satisfies the conditions Aqp > 0, f

< 0, and A(pf = 0.

Now, by invoking isotropy one can conclude that there exists an equivalent

function f = f(P1,P2,P3,2) such that


f A), (4.44)
S

The function can then be used in (4.43), together with (4.35), to obtain








3 (
betr = )2 exp2Acpf n(A). (4.45)
A=1-

Comparing spectral decomposition of b',t in (4.42) and (4.45), one can conclude that


(A) = tr(A);exp 2A A= 1,2,3. (4.46)


(4.46) states that the principle directions n(A) coincide with the trial principle directions

nA), and that the plastic relaxation equation takes place along the fixed axis defined by

the trial elastic state.

Finally, an additive form of the plastic relaxation equation is obtained by taking

the natural logarithms of both sides of (4.46)2. The result reads

et =e,tr -A A =1,2,3. (4.39)
't AA,t MA

(4.39) represents a linear return mapping algorithm in the strain space defined by the

elastic logarithmic principal stretches. In Kirchhoff effective stress space, a linear return

mapping algorithm similar to that presented in [63] can be derived if one assumes an

elasticity operator aAB from the equation

3
PA = aABeO A = 1,2,3. (4.48)
B=l

The result reads

3
PA = p Ap aAB A = 12,3. (4.49)
B=I









Similarity in form between the standard linear return maps of the infinitesimal theory and

(4.48), (4.49) allows the algorithms for the infinitesimal theory to be preserved and

exploited for finite deformation analysis, with the added simplification that calculations

now takes place in the fixed principal stretch directions.


4.6 Constitutive Law for Fluid Flow

Similar to solid phase one needs to describe appropriate constitutive law for fluid

phase. In this study flow is assumed laminar and generalized Darcy's law is employed to

describe the constitutive relation between relative volumetric flow v of(4.7) and fluid

potential n. Linear constitutive equation is given as

S= -k grad n, (4.50)

where k is the second-order permeability tensor and n1 is the fluid potential, defined in

(4.14). The negative sign in (4.50) implies that fluid always flows in the direction of

decreasing potential. Permeability k is an important soil parameter which depends on

other material properties such as: particle size, void ratio, composition, fabric, degree of

saturation [64]. For most practical purposes, k is assumed to be symmetric, positive-

definite.

For incompressible flow the potential H can be decomposed as


n = n +ne; He = t
Jpwg


(4.51)


If and IT represent pressure and elevation counterparts, respectively. g is the gravity

acceleration constant, 0 is Kirchhoff pore pressure as defined in (3.56). Taking spatial

gradient of (4.51), and using (A.3), one obtains






50


grad I = grad 0 + g (4.52)
Jpwg g

If is measured in the direction of gravity, g/g takes a convenient form of {0, 1, 0)T in

Cartesian space e R3. Thus the variational equation (4.18) for the volume conservation

may be written as

H(,l,w)=v J8WJdV+LBgrad-k- .(Pgrad9- +JgdV-J lQdA. (4.53)
P Pwg










CHAPTER 5
HYPERELASTIC-PLASTIC MCC MODEL


5.1 Introduction

Elasto-plastic models based on critical state formulations have been successful in

describing many of the most important features of the mechanical behavior of soils such

as hardening, softening and pressure sensitivity. The modified Cam-Clay (MCC)

plasticity model of critical state soil mechanics [28] is one of the most widely used

plasticity models because of its practicality and simplicity. As a result, MCC model is

adopted in this study to simulate elasto-plastic response of phosphatic waste clay. Two

important modifications are incorporated to the small strain version of MCC to take into

account large deformation effects. These modifications are hyperelasticity and

bilogarithmic compressibility.

Elasticity models are commonly incorporated into elasto-plastic constitutive

models through a hypoelastic formulation. However, extension of a hypoelastic

formulation to the case of nonlinear elastic soil response could result, in some cases, in

conservative models [65]. In case of small strain formulation, the use of non-conservative

elastic models consistent with critical state theory has been justified by the hypothesis of

small deformation [66]. This argument is unacceptable in the large deformation regime

particularly under conditions of cyclic loading where significant energy can either be

extracted or dissipated from certain loading cycles. On the other hand, hyperelastic

materials are those for which a stored energy function exists, and hence, are conservative.

51









Nonlinear hyperelasticity model with a constant elastic shear modulus is used for large

deformation Cam-Clay model in [67]. Though energy conservative, use of constant shear

modulus can be erroneous since experimental evidence suggests that for soil elastic shear

modulus does vary with effective volumetric stress [60, 68]. Consequently, a class of

two-invariant stored energy functions [68] is employed in this study which includes

pressure-dependent as well as constant elastic shear modulus for special case. A variable

elastic shear modulus leads to fully coupled volumetric and deviatoric elastic responses.

Another limitation in small strain formulation is the use of linear variation of the

void ratio (or specific volume) with logarithm of effective volumetric stress to describe

the hardening response of the soil [28]. This assumption can be justified for small

volumetric strain, which does not hold for large deformation regime. In large strain,

linear void ratio logarithm of effective stress variation can result in a physically

meaningless solution such as the prediction of negative specific volume even at

realistically low values of stresses. In this study, this limitation is addressed by

incorporating bilogarithmic compressibility law, i.e., linear relationship between the

logarithm of specific volume and the logarithm of effective volumetric stress as proposed

in [69-71]. Advantages and generality of bilogarithmic compressibility law are discussed

in a following section.



5.2 Hyperelastic Model

The formation of hyperelasticity is based on the existence of a stored energy

function T = I(8c), where Et is the vector of elastic lograrithmic principal stretches. The








effective principal Kirchhoff stress vector 3 can be expressed in terms of'T (see (4.38)).

Substituting (4.38), the elastic moduli a' e R3"3 can be expressed in tensor notation as

e dpi d2 .
a f.(5. 1)
j j

Assuming y(se)= YsP(e,s), one can use the chain rule to expand (4.38) as


Pi -= (5.2:

Sand are the volumetric and deviatoric invariants of respectively.
EV andFs are the volumetric and deviatoric invariants of c, respectively.


se = e .5;
86v=e6o


,+: llel


(5.3)


e =E e l e e3


where5[ 1 1 ]'. Since


-e
V =;
&e


&e J3


where & = e e then (5.2) can be rewritten in the equivalent form

S= p + Jqna = p +s,
where


p =-=-P-8;


-TI;


(5.4)


(5.5)


(5.6)


s = p- p.


p and q are the mean normal stress and the deviatoric invariant of 3, respectively. s is the

vector of deviatoric principal Kirchhoff stresses.


)


)








The elastic moduli tensor can be obtained by differentiating the stress equation

(5.6) with respect to the corresponding strain components. In order to do that, one would

need De = VVY e R2x2, Hessian matrix of Y.


SDe Dfe 2/&ee a2/e e1
De=D D%[821F a12iva c s] (5.7)
LDei De2. a /ee a2 /sese

The first time-variation of stress invariants now takes the form


j= De }v (5.8)


Note that D" is symmetric provided that the function W exists. If De2 0, then the

volumetric and deviatoric elastic responses couple, that is, an imposed volumetric strain

produces a shearing stress response, and vice versa. The following section investigates

the coupled elastic responses within the context of stored energy function developed

specifically for cohesive soils.

Consider a class of stored energy function of the form [68]


SPo p --Evo 3 2 (5.9)


where evo = elastic volumetric strain corresponding to a mean normal stress of po; K =

elastic compressibility index; and t = elastic shear modulus defined by the expression


p= go +apo


(5.10)









p contains a constant term po and a term that varies with ev through a constant

coefficient a. If a = 0 and po > 0, then the elasticity model is defined by a variable elastic

bulk modulus and a constant elastic shear modulus.

The following elastic constitutive equations can be derived from (5.7) and (5.9):


p = oexp v vo ; q= 3pS, (5.11)



where 3=1+3al s 2/2c.

re e e
De =K=- Poexp v vo (5.12a)
K K K
X 1

De2/3=P= 0 o+ ( P=o +aPoexp v evo ; (5.12b)


D e e e3 3o -g vo
D2 =D O s -- J p- exp v -o3 (5.12c)


An important feature of elastic soil behavior, presented in (5.12a), is that the elastic bulk

modulus K is a linear function of p. With a suitable selection of parameters, elastic shear

modulus 4 can be made constant or a linear function ofp (see (5.12b)). Since the

coupling terms D2 and D1 can be nonzero for a 0, the elastic shear and volumetric

responses are coupled for a general loading path. In extreme case, when po = 0 and

es = 42K/3a; det(D') = 0, i.e., D' becomes singular. This situation arises when the

stress ratio q/p reaches its maximum value of 43ai /2 [68].









5.3 Plasticity Model

The essential ingredients of a plasticity model are a yield function, a flow rule and

a hardening law. Two-invariant yield function of MCC model [28] is given by the

ellipsoid


2
f= f(p,q, pc)= +P(P-Pc)=0.
M2


(5.13)


Here f is defined in the space of principal Kirchhoff stresses, P. Invariants p, q are given

in (5.6), Kirchhoffpreconsolidation pressure pc is a plastic state variable that describes

the size off. M is the constant slope of the critical state cone in the p-q plane.


Yield
Surface


Figure 5.1 Yield surface of the MCC Model in p-q plane



Hardening law and flow rule of MCC, appropriate for large strain formulation, are

presented in the following.









5.4 Hardening Law

In case of small strain, the growth ofpc is conventionally defined by a linear

relationship between the void ratio e and the logarithm of pc, or, equivalently, by a linear

variation of the specific volume u = 1 + e with the logarithm ofpc for virgin loading (see

Fig. 5.2a). Corresponding hardening law takes the form

=6 c (5.14)
uo Pc

where uo is the reference initial specific volume at a preconsolidation pressure pco, and 1

is a constant compressibility index of the soil. Upon integration, the hardening law (5.14)

defines the following relationship between the specific volume u and the logarithm of p


uo IPco
=1 !I (5.15)


The limitations of this hardening law are generally well recognized, and include

among others, that a negative void ratio can result even at realistically low values of

preconsolidation pressure, and that the linear relationship is valid only over a narrow

range of values of the effective volumetric stress.

An alternative hardening law for finite volume changes, which appears to have

been first proposed by Hashiguchi and Ueno [69], and later studied more extensively by

Butterfield [70] and Hashiguchi [71], is of the form

6= -1Pc, (5.16)
o PC

where X is the appropriate compressibility soil index in the large deformation regime for

virgin loading. A simple integration of (5.16) yields the relationship










In = -.iPn (Pc (5.17)
o( IPcoJ

which indicates a linear variation of In u with In pC (see Fig. 5.2b).



u In u

UO ---t UO --

I I
S--- ----- ------------- -------------------





po p, In pC P p In pC
pCo PC Pao Pc


Figure 5.2a Unilogarithmic compressibility Figure 5.2b Bilogarithmic compressibility
law law



(5.17) can be also be written in the form


uo Pc


which implies that u -, 0 as pc oo. Since practically one cannot have u < 1 (or e < 0),

the bilogarithmic compressibility law is not without limitation either. However,

Butterfield [70] shows from compression test data on natural soils, specifically soft clays,

that this law is more accurate than the unilogarithmic compressibility equation over a

wide range of values of effective volumetric stress. Furthermore, the value of pc below

which u 2 1 (or e 2 0) is higher with the bilogarithmic compressibility law (see Fig. 5.3).









A simple inspection shows that in the limit of small strains, the natural volumetric

strain In(u/uo) = In(1 Au/uo), where Au = uo u, coincides with the natural volumetric

strain Au/uo of the infinitesimal theory. Thus, the bilogarithmic hardening law

approaches the unilogarithmic law in the limit of small volumetric strains. However,

large deformation analysis requires the use of natural, and not nominal, strains, and so

(5.16) is more robust since it is useful both for small and large deformation analyses. In

light of these desirable features, bilogarithmic law is adopted in the proposed model.



u

DO--




Eqn. (5.16)
1.0 -------- --------- -
Eqn. (5.14)


Poo Po

Figure 5.3 Limit of validity: comparison between unilogarithmic and
bilogarithmic compressibility laws



In order to develop a hardening law appropriate for large strain plasticity model,

one needs to exploit the properties of deformation gradient F. The product decomposition

of the F given by (4.27) produces the identity


J = det(F) = JeJp,


(5.19)









where Jr = det(FV) and J = det(FP). In the space of principal logarithmic stretch,

3 3 3
ev = In J=n(ai2X3)= PSA; E =InJe = se and P = InJP = es Natural
A=I A=1 A=I

logarithm of (5.19) then yields

In J=Ine +lnP =>ev =ev +es. (5.20)

In other words, the product decomposition ofF is equivalent to the additive

decomposition of the natural strains (see Section A.3). The rate form of (5.20) is

1 jie JPe
==-> = EV +ey. (5.21)
J J J

J/J = 6/u since I = u/ uo. u is the specific volume with a reference value uo in the

undeformed configuration. Thus, the hardening law of(5.16) can also be written as


S = ev + E v- -. (5.22)
J u Pc

Setting se = 0, q = 0 and p = pc in (5.11) yields the following expressions for virgin

isotropic loading:

e e
Pc = Poexp Ev v ; (5.23a)


e C
K I V Pc Ve
tc *e e -. (5.23b)


Substituting (5.23b) in (5.22) and simplifying, one can obtain the following hardening

law expressed in terms of plastic component of the natural volumetric strain:








Pc= Of; = (5.24)
Pc ,-K

Integration of (5.24)1 produces

Pc = Pc,n exp[O(p, -eSn)] = Pcn exp[o('tr -e)], (5.25)


where pc,. and ep,n are the preconsolidation pressure and the plastic natural volumetric

strain at time tn, respectively. Thus, the hardening law given by (5.16) offers a further

computational advantage in that the evolution equation pc can now be integrated exactly

over a finite load increment.


5.5 Flow Rule

Under the hypothesis of associative flow behavior, the integrated flow rule at any

time t.n+ in the space of logarithmic principal stretches takes the form (cf (4.47))

,e =.e,tr _A (5.26)
ap,

For simplicity of notations, in (5.26) and in the following subscript (n + 1) is omitted; it is

assumed that the unsubscripted variables are all reckoned with respect to time station tn.j.

Volumetric and deviatoric components of (5.26) are as follows:

,e =evtr A ; ee = eetr -A, Of s (5.27)
8p C2 v t 1 taq

In terms of unit vectors ^, tn (5.27)2 takes the form

e etr tr 8,f (
el =B n -A -n, (5.28)
Bq








where es =e -/e 1 etr = tr A tr =etr_ /tr = ee/lee. Two useful

identities, derived from the following assumptions, are: (i) from the assumption of

associative flow rule, n = ee/i ee = s/s.ll and (ii) from the assumption of convexity of

the yield function, f = itr. Exploiting these identities (5.27) can be rewritten as

se= etr Ae; e e= se tr A( (5.29)

subject to the conditions

f(p,q,pc)0; A> 20; Agf(p, qpc)=0. (5.30)
tr tr e, tr e etr
In the elastic regime, f(p q ,Pcn) < 0, A( = 0; e = sv ,e = e Plastic regime

is realized for the conditions: f(pt qt, Pc,n) > 0, Aq > 0. ptr and qtr are the predictor

values defined from (5.6) as


&e,tr q &e,tr

In the elasto-plastic regime, (5.29) and (5.30) can be viewed as a system of simultaneous

nonlinear equations in the elastic strain invariants and the consistency parameter A(p,

represented by residual vector r and vector of unknowns y, as follows:

ev -Evtr +Aep f/sy
r= l-s'r+A f/aq ,; y={ } (5.32)
f Aq

In order to solve this system iteratively, one can employ Newton's method over the

following loop:









yk+l yk (klrk; k a rk
Sk+ k I k; Ak- (5.33)
y k

k is the local iteration counter. AER3"3 is the consistent tangent operator. A closed form

expression of A can be written in a compact form with the aid of the following matrices:

H=H, H1 H12 _Ja2/fapap 82f/apaq]
H21 H22 2f/la p a2//aqaq (5.34)
G=HDe.

H = VVflpc e R2x2 is the Hessian matrix of yield function f with pc held constant. The

A matrix then takes the form

S l+Aq(Gl+Kp a2/lapapc) AqpG12 Of/ap
A= AGG21 l+AqpG22 Qf/ q (5.35)
De la / ap + D la /aq + Kpa /apc D e2a / ap+ D 2af / aq 0
DiIf/8p+De9f/9q+KpOf/Opc D18f/p+D4 f/8q 0

where Kp = Opc /&ev = -Opcis the plastic hardening modulus. Using the yield function

(5.13), one can have: af/ap = 2p- pc, f /aq = 2q/M2, a/f pc = -P,

2f /pPpc = -1; H11 =2, Ha = 2/M2, H12 = H21 = ; G1 = 2Dfi, G2 = 2DDe2

G21 = 2De, /M2, G2 = 2De /M2.



5.6 Consistent Tangent Moduli
Material tangent stiffness matrix aeR3x3, defined in the spaces of and e, is
expressed as

(5.36)
ai &j = ,tr
J








For a fully elastic response with volumetric and deviatoric coupling, the matrix a takes
the form

a= D- q 55+ D(e (853)ii +-8)+ 2q (I- i)+D @i, (5.37)
9e J3 3e 3

where I is a 3x3 identity matrix, ii is a 3x1 vector obtained from the relations

S=ee/leel =eetr /etr = /s/ll.

In case of isotropic, linear elasticity free energy function Y takes the form

Y=1 X[, + + Fe + eS 12)2+ )2 )2
2 (5.38)



where K and p are constant elastic bulk and shear moduli, respectively. K = X + 2/3p, X

is a Lam6's constant. Since p = po > 0 (see (5.10)), a = 0; elastic shear and volumetric

responses uncouple, i.e., De=De = 0 (see (5.12c)). Now substituting De, = K, De=

3p (see (5.12a), (5.12b)), tangential elastic moduli of (5.37) degenerates to the familiar

expression for linear elasticity as follows:

ae =K6 6+2 I- 6@8). (5.39)

In the elastoplastic regime, the tangential moduli matrix a' can be expressed

using strain derivative of (5.5) as

aep= a 8 p J2 q+J2 0 (5.40)
&e,tr &e,tr I3 etr (e5tr


where









A a e/etr eetwr/e e3 1
&e,tr &etr &e, tr ee, tr (3


(5.41)


Substituting (5.41) in (5.40), and using the elements of the matrix D* to enforce the chain

rule, one can have


Df


3 0 e t r
3e Dt2e-+ D t-
J3 I atr &eatr


(5.42)


+ rI-2q -
3,etr 3


Strain derivatives of the invariants ev and Ee are obtained from (5.29) as


v (e,tr
&e,tr &e,tr v


-Ap4 P
ap )'


- e,tr A .(5.43)
&e,tr &e,tr s cJ)


In the expansion of (5.43), one will need the following strain derivative of pc


Pc,
e, tr


= K +Ktr
e, tr + p


(5.44)


where Kp = 8pc/e and K = apc/ r The expansion of (5.43)1 and (5.43)2 takes

the following forms, respectively


e a f CA (
+bz12 S -S
&e,tr ap -e,tr


&
+b22zz -
&e,tr


Of BeAq
aq ae~tr


Ctr
&e,tr


b2 e,
&e,tr


where


(5.45a)


(5.45b)









bl1 =I+A{ Gi +Kp ;


b21=A qiG21 + Kp ;2f)


c= -ApKtr a2f


Solving (5.45a) and (5.45b) simultaneously yields


E b22C1 5- bl2n
&e,tr det(b) 22

ase! -1 bii c+ b
atr det(b)21c+ 3b11i
&e,tr det(b) 3


b12 = AqG12;


b22 = l + AqG22;


- b22



aq


a f CIA P
-bg I
-b12 q e, tr



b21 ~Ie,tr'


where det(b) = blb22 b21b2.

The strain-gradient of Ap is obtained from the overall consistency condition


_ f Op
op &etr
p g e,tr -


Of Oq +
aq etr


Of Pc
fPc at = 0.
opc &e,tr


Since p, q and pc are functions of the strain invariants, one can expand (5.48) further by

chain rule, and then use (5.47a) and (5.47b) to obtain the following result


Aq = a8 + a2ii,
e, tr


al =[dlb22c-d2b21c+det(b)Ktr ,
e O p


(5.49)


a2 = 1 (d2b I-dib12 ,
e


-b12f)


+d2 b, a-


d e =D +Del +Kp af
d, =D1j +D21 +K t
ap Oq PaP)


d e af
d2 D12-+
4p


(5.46)


(5.47a)



(5.47b)


af
aEe,tr


(5.48)


where


b2 1'


(5.50)


De f
D22 `
Oq


e= dl b22








Plugging (5.49) in (5.47) yields


a&
e = DP5
e, tr


fip2 ;
3 12


8e, =_D 3 p
n'tr 21 3 22


(5.51)


Dp e R2"2, consists of coefficients of base vectors 5 and f in (5.47a) and (5.47b), is

defined as follows:


DP = )[b22 c-
det(b) I


al pL+bl2al ,


D det(b) bl2( 1+ a2 b22 2 ];
12 2-de2t(b) I J2 N

DP- de=(b) bllal '" b21 ac l I
D2P 21 C- 81 a 3
D det(b) [b1 8q 2a I 2


The strain gradients of the stress invariants then take the form
The strain gradients of the stress invariants then take the form


ep = Dep5+ D epi;
&-,tr IpI8+ 3 n


=q DeP
&etr 21


where D" e R2"2 is defined as


Dep = DeDP.


Substituting (5.53) in (5.40) then yields consistent elasto-plastic tangent moduli


aep =D ep -s
s
+ 2q (I-i@i)+
3s


+ jDit2


(5.55)


2 D~e@ii.
322


For elastic loading, DP= I, D' = D', and so (5.55) degenerates to (5.37). Thus (5.55)

represents a generalized expression for both elastic and plastic loading. For elasto-plastic


(5.52)


+ D -ep
3 22 n,


(5.53)


(5.54)






68


loading DP # I, and so a" loses its major symmetry due to the fact that Dep D .

Volumetric and deviatoric responses are coupled for elasto-plastic loading even if D' is

diagonal (in case of constant elastic shear modulus, i.e., po > 0, a = 0), since the matrix

D'P is generally full due to plastic volumetric and deviatoric coupling inherent in the

Cam-Clay model.









CHAPTER 6
HYPERELASTIC-VISCOPLASTIC MCC MODEL


6.1 Introduction

Hyperelastic-plastic MCC model is further extended for viscoplasticity to model

time-dependent secondary compression of phosphatic waste clay. Elasticity response is

based on the stored energy function [68] of(5.9). Consequently, nonlinearity of elastic

moduli and coupling of volumetric and deviatoric elastic responses follow the same

constitutive relations as presented in Section 5.2. Yield function of MCC (see (5.13)) is

coupled with a time rate flow rule to simulate viscid response.

Clay is a strain hardening, rate sensitive material that has remarkable

characteristics such as rate sensitivity of strength, secondary compression, creep and

stress relaxation. Various elasto-viscoplastic constitutive models have been proposed to

describe the theological behavior of clay. Most elasto-viscoplastic constitutive models

can be classified as overstress models or non-stationary flow surface models. Overstress

elasto-viscoplastic constitutive model was first introduced by Perzyna [72]. The

Zienkiewicz et al. model [73], Adachi and Oka model [74], Dafalias model [75], Katona

model [76], Baladi-Rohani model [77] belong to this category. The key assumption in

these models is that viscous effects become pronounced only after the material

undergoes yielding, and that viscous effects are not essential in the elastic domain.

Overstress model is an outgrowth of classical plasticity where viscous response is

introduced by a time rate flow rule with a plastic yield function. As opposed to

69








overstress model, in non-stationary flow surface model yield condition of a material

changes with time as plastic straining occurs. Olszak and Perzyna [78] initiated this

concept by introducing the time dependent yield condition. Later Sekiguchi [79],

Dragon and Mroz [80], Nova [81], Matsui and Abe [82], among others, adopted this

concept. Viscoplastic rate equations of the non-stationary flow surface model are

characterized by the stress rate terms.

Overstress viscoplasticity model is adopted in this study to simulate secondary

compression behavior of phosphatic waste clay. Motivations for selecting overstress

model are: (1) the incorporation of MCC yield function is straightforward; (2) the

generality of time-rate flow rule offers the capability of simulating time-dependent

material behavior over a wide range of loading; (3) the formulation is amenable to finite

element implementation.



6.2 Flow Rule

In viscoplasticity formulation, additive decomposition of t takes the following

form

t = e +vp, (6.1)

where e is the vector of principal logarithmic stretches. e and e" are the elastic and

viscoplastic components of s, respectively.

For an associative flow Jvp is given by the relation

tvp = y(f)6 (6.2)
MP








where 0 is the vector of principal Kirchhoff stresses, y is a material property called the

fluidity parameter (units of inverse time) that establishes the relative rate of viscoplastic

straining, p(f) is a scale function dimensionlesss) of plastic yield function, f.

,(f) = f> (6.3)
0, f0

p(f) is called viscous flow function. Two commonly used forms of (p) are:


<(/) =- ; (6.4a)

( N
p(f)=exp (- -1. (6.4b)

N is an exponent; fo is a normalizing constant with the same unit as f so that p(f) is

dimensionless. Although more elaborate functional forms of ((f) may be established,

the forms given by (6.4) appear to suffice for many geologic materials [73].

(6.2) can be written in incremental form as

AvP =YAt(p(f) =Ay(f) (6.5)
op op-

where Ay = yAt. Substituting (6.5) in additive decomposition of natural strains, i.e.,

a = se +6P = e,tr + p one can write the flow rule as

se = e,tr Ayq((f) (6.6)
6A )CID (6.6)

Volumetric and deviatoric components of (6.6) are as follows:

e =setr Aa-f; ee = eetr- A (/ (6.7)
vop 2 pq |s||








In terms of unit vectors i, tr (6.7)2 takes the form

aef= e,trftr AY ) fi, (6.8)


where ae = 2-II3e e tr =tr tr= etr/= tr ,f =ee/e Since ii= tr

(see Section 5.5), (6.8) can be expressed in terms of scalar coefficients. Consequently,

(6.7) can be rewritten as

t =s -A () ; e,tr -Aycpqf)L. (6.9)


In the elastic regime, f(ptr qtrPc,n), <0 = = v ,s = s etr

Viscoplastic regime is realized for the conditions: f(ptr,qt,pcn ) > 0, W(f) > 0.

ptr and qtr are the predictor values as defined in (5.31). In the viscoplastic regime, (6.9)

can be viewed as a system of simultaneous nonlinear equations in the elastic strain

invariants, represented by residual vector r and vector of unknowns y, as follows:

r= I -VSt + AY p(f) af /ap s (6.10)
ae -str + Ayq~f)f/f ej

One can employ Newton's method (cf. (5.33)) to solve this system iteratively while the

tangent operator AeR2x2 of (5.33) now takes the form

'An At Brt /& 1r, /&e
.A[A A12 1 [1r (6.1 1)
LA21 A22J &r2/&e &2DB
La~ v s].,/a


Individual components of A are derived as








All =+Ay qf GI +Kp 0 + '(f) Dfl
Ip~c ap OP1a


+D


Al2 = f ')De afL
[(P,(f)Lf ( OIp


A21 = Ay[q'(f)f (De
2^A T D I


A22 =1+ AYq)'(f)ID f +
B Op


De1 + P(f)G2


(6.12)


+D21 +
aq


D + q(f))G] .
d qL


Kp = dpc / ae = -Opc (see (5.25)) is the plastic hardening modulus,

q'(f) = a&p(f)/8f. Matrices D', G are defined in Sections 5.2 and 5.5, respectively.


6.3 Consistent Tangent Moduli
Consistent tangent moduli matrix in elasto-viscoplastic regime a'P ER3x3, defined
in the spaces of J and 8, is expressed as


ep
a.
IJ


Sapi ap0
&j etr


Following the developments of Section 5.6, one can obtain an expression of a',
identical to (5.42), i.e.,


aep =8


(9@ DJ0- e +De


(6.13)


+ 3-tr2q (I-5@8-n .
3 e, tra 3

Strain derivatives of the invariants es and ss are obtained from (6.9) as


e,tr &e,tr (ser


aE8 e,8 tr AT r)J).
ae,tr & e,tr '-aqJ


(6.14)


-f )


Kp +p(/f)G2, ;
Pc )









In order to expand (6.14) to the lowest order, one will need the following strain

derivatives:


_a = al+a2 e+a3 s
ae,tr e, tr ae,tr


= 846 + 5-- +G12 s+
apce, tr e, tr e,tr'

a2f ce ace
=r G21 e+ G22
Lqc&e, tr etr e, tr


(6.15a)



(6.15b)



(6.15c)


Coefficients of (6.15a) and (6.15b) are given as


tr 1
al = Kp ;.
1 P c P

4 Ktr c2
Sapapc


e 81
a2 =D11 +
op


a5 = G +Kp


e eI
D21 +
-+


K Of
CP c


e D a
a3 = D12
9p


921
OpOpc


= -Op and Ktr = 8pe etr = Opc (see (5.25)).


a&e
e,

ev
-e,tr


where


bl = 1+a2Aycp'(j)- + Aycp)G21;
pb21
b2l = a2AY(P'(J)-+ AY(Pf)021;
Nq


ci =- -aIAYp'(f) ;1
op


L+b12 s -I1,
tr a+e,tr


b22e, -28 + -c 3i,
Se,tr 3



b12 =a3ATy9P(f)- +Ay()G22;
ap
b22 = 1 + a3Ay'(f)- + Ayp(f)G22 ;
oq


C2 = aAY(p()-.
aq


where Kp


+De af.
D22 a'


(6.16)


(6.17a)



(6.17b)


(6.18)


= apc ev








Solving (6.17a) and (6.17b) simultaneously yields


e DP5 + DPn; &DP +
e,tr 1 3 e2 e,tr 3 22n

Coefficients ofDP' R2"2 are given as


DP e = b) (b22C1 +b12c2);

1de= ( (-b21C1 b1l2);
det(b)


12det(b)1

D = d1 b1l,
Sdet(b)


where det(b) = blb22 b2lbl2.

Substituting (6.19) in the expression of a' (cf. (5.42)) then yields consistent

elasto-viscoplastic tangent moduli as follows:


aeP =


@ii+ iDep@8
Ofi+ J3 12


(6.21)


+ -4(I-n@i)+2DePi i.
3es 3


Dep = DeDP e R2x2. Note that (6.21) and (5.55) have identical expressions. (6.21)

represents a generalized expression for both elastic and plastic loading. For elastic

loading, D' = D' since D = I and so (6.21) degenerates to an expression of a' identical

to (5.37).


(6.19)






(6.20)










CHAPTER 7
LINEARIZATION


7.1 Preliminaries

Some useful formulas are summarized below. These will be helpful for

linearization of strong and weak forms of coupled equations (see Chapter 4).

The first of these formulas is the Piola transformation, introduced first in Section

4.2. Let y E Rmd be a vector field in spatial configuration h(B). Then, the Piola

transform ofy in reference configuration B is

Y = JF-- y, (7.1)

provided that motion < is regular in B. The following equation holds for y, Y.

DIV Y = Jdivy. (7.2)

Proof of (7.2) is given in Section A.7. This identity may be extended to cases where Y

and y are vectors derived from tensors of order greater than or equal to two by fixing all

but one of the tensor's legs (for example, fixing one leg of the Kirchhoff stress tensor r

produces a vector of Kirchhoff stresses).

Following are linearization of some basic terms, one would need for subsequent

development. Let 8u be the variation of the displacement field; then the linearization of

F and F' at any configuration 4(B) are given, respectively, by

LF= F + grad8u F= F + GRADSu; (7.3a)

LF"- = F"' F-1 grad 8u = F"1 -F"- GRAD 6u- F1. (7.3b)

76






Derivations of(7.3a) and (7.3b) are given in Section A.8. Linearization of the Jacobian

and the rate of the Jacobian at spatial configuration t(B) are given, respectively, by

LJ = J + Jdiv(Su); (7.4a)

LI= + J[div(8v) -gradv: gradt(8u)+div(u)div v]. (7.4b)

See Section A.9 for derivation of (7.4a) and (7.4b).

Linearization of the reference saturated mass density po = Jp at the spatial

configuration t(B) is

Lpo = Po +PwJ div(8u). (7.5)

Proof of (7.5) is given in Section A. 10. Note that po is not constant since, as pointed out

previously (see Section 3.5.2), the total mass of the soil-water mixture in B is not

necessarily conserved in Wt(B). The variation of po reflects the amount of fluid that enters

into or escapes from the soil matrix due to the variation of the Jacobian.



7.2 Linearization of Strong Form

The results discussed in the previous section can be applied for the lineraization

of strong form or the field equations of linear momentum and mass conservation. Since

the field equations are mixed formulation involving finite deformation u and Kirchhoff

pore pressure 8, linearizations are derived consistent with the imposed infinitesimal

variations 8u and 80.








7.2.1 Equation of Equilibrium

Let E = DIVP + pog be the linear momentum equation (see (4.3)). Substituting

(3.55) and (3.56), E can be rewritten as

E = DI +OF-t)+ Jpg. (7.6)

Taking the variation yields

E = DIV(SP + F-t + 0 F-)+ 5(Jp)g
BE = DI (7.7)
=DIV(A :F +0F-t +SOF-t)+B(Jp)g.

A is the first tangential elasticity tensor of order four. A has the structure Aijd = P1ij/OFu.

One can write from (7.5)

6(Jp)g = pwJdiv(Su)g. (7.8)

Now, using (7.1) and (7.2) one can express

Jdiv(Su) = DIV(SU) = DIV(F-1 Su (7.9)

where BU is the Piola transform of Bu. Substituting (7.9) in (7.8) results

8(Jp)g = pwJdiv(Su)g= pwDIV(JF-1 -5u)g. (7.10)

Substituting SF, SF' from (7.3), 5(Jp) from (7.10) in the expression of BE of(7.6),

linearization of E may be written as

LE=E+DIV(A: GRAD u)-DIV(OF-t -GRADtSu F-t)
+DIV(68F-t)+pwDIV(JF- u). (7.11)

A crucial step in the linearization of the linear momentum balance equation is the

evaluation of the tangential elasticity tensors of the solid matrix. Four of them are

introduced in this section: the tensors A, D, a and d. Each of these tensors can be derived









directly from others. For example, A can be obtained from D using the following

expression

A= 2F.DFt +S 1; Aijk = 2FimFknDmjnl +SjlSik. (7.12)

See Section A. 11 for derivation of(7.12). D = 8S/DC (i.e., Dmjn = aSmj/oCi) is the

second tangential elasticity tensor of order four. Constitutive relation (4.36) and pull

backs of the Kirchhoff effective stress tensor r yield the following expressions of the

second Piola-Kirchhoff effective stress tensor S as

3
S=F-' -F-t = XPAM(A); M(A) =F-1.m(A) F-t. (7.13)
A=I

PA'S are the principal Kirchhoff effective stresses and m(A)'s are the dyads formed by

juxtaposing the principal directions of the elastic stretches, as given explicitly by (4.35)2.

Using the chain rule, from (7.13)i one can obtain the following expression for the tensor

D as


D=-=S 1 apM A) M(B)+ pA ~ (7.14)
A=IB=1 B A=1

since &s/OC = M13)/2 (see Section A. 14).

By the symmetry of both S and C, and by the axiom of material frame

indifference, the tensor D possesses both the major and minor symmetries such as Dij =

Djiu = Diyn Diu = Duj. Spatial tangential elasticity tensors a and d may be obtained

from the push-forwards of A and D as

aijkl = FjaFlbAiakb; (7.15a)

dijd = 2PaFjbFkcFldDmw. (7.15b)








a and d have the symmetries [83]: aijk = auij, diju = djua = diju = duij. A push-forward of
all the indices of D gives the following expression for the spatial tensor d

33 3
d = i(A) @m(B) + 2pAd(A) (7.16)
A=1B=1 A=1
Here d(A) is the push forward of M(A)/d as

aM(A)
d A = FiaFjbFkcFld ab (7.17)


For the general case of left Cauchy-Green tensor b having distinct eigenvalues )X, 2

and 1 d(A) takes the form

d(A) =-- [Ib -b@b+I3A2(l1-I)
DA
+- A [x (b Om(A) +m(A) Ob)- IDXAm(A) m(A) (7.18)
DA 2
13 A (l "m( +m (A) ]
DA

where (Ib)ijkl = bikbjl +blbjk],I3 =det(C), Iijkl =6kikil +6ilijk

DA = ( A ) ), D'h A = 8A ZIXA 213A3. {A,B,C} denotes an even

permutation of the indices {1,2,3). This expression is strictly valid only if the
eigenvalues are different since DA = 0 otherwise. Although it is possible to derive a
similar result for the case of repeated eigenvalues [19], from a computational standpoint
it is preferable to reduce this situation to the general case of distinct eigenvalues by
introducing a perturbation of the repeated roots. For example, in case of repeated roots, a








small perturbation of 10 x repeated eigenvalues is used in PlasFEM to obtain a general

case of distinct eigenvalues.

The spatial counterpart of (7.11) may be derived directly from Piola

transformation. Then, the linearization of E in the spatial picture takes the form

LE=E+div(a:grad6u)-div(Ggradt6u)+grad(80)+Jpwdiv(5u)g. (7.19)

An equivalent form to (7.19), using the spatial tangential elasticity tensor d, is

LE= E+div[(d + T' ):grad6u]-div (gradt8u)+grad(68)+Jpwdiv(5u)g, (7.20)

since a = d + T l 1(see Section A. 12). Here (r S 1)ljld = Tjl8ik represents the initial

stress contribution to the spatial stiffness. Comparing term by term, equivalence between

(7.11) and (7.19) can be established. Exploiting the Piola identity of(7.2), one can show

that DIV(A:GRAD Bu) = div(a:grad 6u) since A:GRAD Su is the Piola transform of JF

a:grad 8u and grad J = 0 (see Section A. 1). Similarly, second divergence terms can be

shown equivalent noticing that 0 F1. GRAD' u F'= 0 grad' u F' is Piola transform

of T'(0 grad' 5u). Expansion of the third divergence term of (7.11) yields DIV(50 F') =

GRAD 58 F' + 60 DIV(F) = grad 80 since DIV(F4) = 0. Equivalence between the last

terms can be derived from (7.9).


7.2.2 Equation of Flow Continuity

Let M = DIVV + DIVV be the equation of flow continuity for a saturated soil-

water mixture in material reference (see (4.9)), where V and V are the Piola transform of








v and V, respectively (see Section 4.2). Then, the linearization ofM at any configuration

(B) is

LM = M +DIVSV +DIVSV. (7.21)

One needs expressions of BV, 8V for linearization of M. First consider variation of V,

6V = JF-1 v)= 5JF-1 v + JF-1 v+JF-1 *-v. (7.22)

From (7.4a), (7.9) and (7.3b), one can express

8JF-i-v=DIV(JF-1 8u)-1 (7.23)
( '(7.23)
JF-1 v =JF-1 grad u- v = -JF-1.GRAD8u F-1 v.

Substituting (7.23), one can rewrite (7.22) as

V = 5(F-1 v)= DIVJF-I1 .u v-JF1 GRADBu F- + JF-1 8v. (7.24)

From (4.50) and (4.52), generalized Darcy's law can be written as

Jg= -k P +i (7.25)
1Jgpw gJ
Piola transform of V yields

V= JF-.= -K- 0 A +JFt. (7.26)
I gPw 9)

where K = F-1 k F-t is the pull-back permeability tensor. K and k are assumed

symmetric in following derivations. Using chain rule, one can expand 8V as

(GRADG GRADES
5V=-BK- GRDO +JFt -K. GRAD Ft).J. (7.27)
S gPw g gPw
Substituting 6J from (7.23)1 and 8F from (7.3a), one can write








6JFt)= DIV(JF-1 .uFt + JGRADt8u. (7.28)

Substituting (7.28) and rearranging, (7.27) can be rewritten as

8V = -K -GRAD5 + [DIJF-1 -8u)Ft + JGRADt5u 1
(7.29)
-SK- (GRAD0 +J.IF
G 8Pw t

where

6K=-F- .{2Sym(GRADu -K-Ft)-(I+eo)DIV(JF-' .u)-ak}.F-t.(7.30)

If permeability tensor k is assumed to be varying with the deformation or void

ratio of soil skeleton variation of k can be expressed as

8k = (+ eo) J div(u)& = (1+eo)DIV (JF' u) (7.31)
&e ae
where eo is void ratio of the soil-water mixture in reference configuration. See Section
A.17 for derivations of (7.30) and (7.31).

Using (4.8) and (4.17), M can be written in spatial description as

M= Jdiv v + Jdiv = +Jdiv (7.32)

Corresponding lineraization takes the form

LM= M +(J +Jdiv))
= M + +8Jdiv V+ J (div V).

Substituting (A.31), JS(div V) can be expressed as

J (divV) = Jdiv(58)- grad(JV): gradtsu, (7.34)








since JgradV = grad(J) by knowing that grad J = 0. Substituting (7.34) in (7.33) then

yields

LM = M + 8 + div[5(J)]- grad(J): gradt8u, (7.35)

since

div [8(Jr)] = div (8J V) + div (J 8v)
= grad J -V + 8J div + gradJ 5 + Jdiv(8-)
= (GRAD 5J F-1). + 8J div V + J div(8v) (7.36)
= 5(GRADJ)- F-' + J div + J div(8v)
= 8J div V + J div(50).

Component terms of (7.35) are given as

8i= J[div(Bv) -grad vgradt (u)+ div(Su)div v; (7.37a)


JV = -kI grade g (7.37b)


8(J) = -k grad(80)- grad grad58u Jdiv(
gp(7.37c)
-(l+eo)Jdiv(8u)k.[grad +j .
S Pw

(7.37a) and (7.37b) are obtained directly from (7.4b) and (7.25), respectively. Derivation


of(7.37c) is given in Section A. 19.








7.3 Linearization of Weak Form

Linearization of weak form G(O, H, 1I) of(4.13) takes the form

LG=G+B gradr:(d+TrE1):grad5u dV+fJ (8div-08gradti:grad8u)dV

-f pwJdiv(8u)l:g dV- Tl-.6t dA.

(7.38)

where 8u, 56 and Bt are the variations of the displacement vector, Kirchhoff pore

pressure, and traction vector, respectively.

First integral term of (7.38) is derived from the variation ofgrad rn: r (see Section

A.20). The variation 6(8 div r1) = 88 div rI +8 65(div 11) produces the second integral term

upon substitution of the identity 8(div 1)=div 8n gradtl : grad8u = -gradtt : grad8u

following (A.3 1) (note that 86n = 0 since I is a vector of arbitrary virtual displacements).

The third integral term emanates from the linearization of po (see (7.5)). The last integral

term is derived from a straight-forward linearization of the traction vector t.

Upon substitution of (3.55) and (3.56), linearization of weak form G(, I, Tq) in

material description (see (4.12)) can be expressed as

LG = G+ f8(GRAD1 T: P-poT g)dV- 5(Ti.t)dA
=G+J 8(GRAD q:P)dV+J 8(GRAD : (OF-t))dV (7.39)
J8(Po1i g)dV- J'8(T. t)dA.

Variation of the first integrand is given as 8(GRAD -T: P) = GRAD 1q: A: GRAD 8u (see

Section A.21).Variation of the second integrand in material description may take the
form








8(GRADTi: (bF-t))= GRADI :F-t -8 GRAD : (F-t GRADtSu F-t)

Upon substitution of (7.10), variation of the third integrand yields

8(poT- g)= PwDIV(JF-1I u)I: g.

Substituting all these identities in (7.39), one can obtain linearization of G in material
description as

LG=G+ GRADTi:A:GRADu dV+ I 8GRADrl:F-tdV
-IB -RAD: -t(.GRADt8u-F-t)dV (7.40)

-f pwDIV(JF-1- 8u gdV- T-8tdA.

(7.39) and (7.40) are equivalent expressions.
Linearization of weak form HAt((,fI, n ) of (4.25) is given as

LHAt =HAt +f Jdiv(Su)dV+PPf gradW- k grad59 dV
B At gpw

-2pPo0 gradJ -Symr k gradtsu) grad0 dV
Js LgPw J

-P33o gradW-[gradu-(div 8u)1]. k .J dV (7.41)
g

+PoO(l+eo) gradv.(div Su)l. -k + j1gr JdV
J B lPw8 8g
-PPo f313j Q dA.

The first integral term of (7.41) results from the variation of J (see (7.4a)). The second,
third, fourth and fifth integral terms result from the variation of grad v- J (see Section

A.22). The last term results from the variation of Q.
Linearization of weak form HAt (, I, w\) of material description (see (4.26)) is


given as




Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EZZBLI29T_5JFZN0 INGEST_TIME 2013-04-08T23:19:33Z PACKAGE AA00013546_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES


xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID ESRBCC5EZ_GZSBXH INGEST_TIME 2013-03-12T14:21:24Z PACKAGE AA00013546_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES