UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Subjects
Notes
Record Information

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 HYPERELASTICPLASTIC 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 HYPERELASTICVISCOPLASTIC 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 Onedimensional 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 Onedimensional 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 Hyperelasticplastic MCC model parameters for large strain simulation of laboratory consolidation tests ................ 125 10.4 Hyperelasticplastic MCC model parameters for small strain simulation of laboratory consolidation tests .............. 125 10.5 Hyperelasticplastic MCC model parameters for simulation of laboratory consolidation test ..................... 126 10.6 Hyperelasticviscoplastic 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 hyperelasticplastic consolidation (pond depth 2.44 m) .................................................. 133 10.11 Material parameters for hyperelasticplastic consolidation (pond depth 4.57 m) ................................................. 134 10.12 Material parameters for hyperelasticplastic 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 hyperelasticviscoplastic consolidation (pond depth 2.44 m) ................................................ 141 10.16 Material parameters for hyperelasticviscoplastic consolidation (pond depth 4.57 m) ................................................ 141 10.17 Material parameters for hyperelasticviscoplastic 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 twophase 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 pq 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 onedimensional consolidation problem .............................. 99 9.3 Onedimensional hyperelastic consolidation: variation of total potential with time ........................................... 100 9.4 Onedimensional hyperelastic consolidation: isochrones of constant Cauchy pore pressure ................................. 102 9.5 Onedimensional 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 Crosssection of clay slime ponds of surcharge area no. 1 (Source: PSI) .................... ............ ........... ... ...... .... 121 10.10 Crosssection of clay slime ponds of surcharge area no. 2 (Source: PSI) ............................... ............ ......... ... .. 122 10.11 FE mesh for oedometer cell ................................................... 123 10.12 Large strain, hyperelasticplastic simulation of laboratory consolidation tests ................................. 123 10.13 Small strain, hyperelasticplastic simulation of laboratory consolidation tests .................................. 124 10.14 Hyperelasticplastic simulation of laboratory consolidation test: large strain versus small strain ............... 126 10.15 Large strain, hyperelasticviscoplastic 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 Hyperelasticplastic consolidation settlement: pond depth 2.44 m ................................................... 135 10.19 Hyperelasticplastic consolidation settlement: pond depth 4.57 m ................................................... 136 10.20 Hyperelasticplastic 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 Hyperelasticviscoplastic consolidation settlement: pond depth 2.44 m .............................. .................... 142 10.24 Hyperelasticviscoplastic consolidation settlement: pond depth 4.57 m .................................................. 143 10.25 Hyperelasticviscoplastic 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 boundaryvalue 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 soilwater mixture as a twophase 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 CamClay (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 displacementbased 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 Nonlinear response ofgeotechnical structures typically results from plastic yielding and large deformation of the soil skeleton. There are many classical geotechnical applications where nonlinear 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 P8 effect. The impact of large deformation and plastic response is most evident in soft clays. It is a wellknown 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 soilwater 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 timedependent: 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 noflow (undrained) and the freeflow (drained) conditions can be analyzed using a singlephase continuum formulation, consideration of a twophase soilwater 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 [110]. 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 twofold 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 rateconstitutive equations [8,9,1214]. 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 longstanding issue in large deformation consolidation analysis. Classical theory of mixtures [2126] is employed in this study to describe coupled response of solid and fluid phases. Accordingly soilwater mixture is viewed as a twophase 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 CamClay (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 nonlinear large deformation consolidation analysis is implemented in a displacementbased 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 solidcontent 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 noninteracting mixture of nonpolar 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 boundaryvalue 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. Hyperelasticplastic MCC (modified CamClay) model, suitable for primary consolidation response, is presented in Chapter 5. Chapter 6 presents the development of a hyperelasticviscoplastic MCC model, used for simulation of timedependent 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 twodimensional 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 multilane expressway around Lakeland, Florida). Numerical simulations are run for cases of hyperelasticplastic (primary) and hyperelasticviscoplastic (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 fluidfilled 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 threedimensional 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 quasistatic 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 fluidsaturated 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 elastoplastic nonlinear constitutive models and obtain reliable solutions of the field displacements and pressures. A general analytical procedure that accounts for nonlinear 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 nonlinear 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 porefluid 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 [2126]. This research is focused on fully saturated porous media (twophase 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 socalled natural formulation. Subsequently, however, these authors appear to favor hypoelastic rate models on the basis that multiplicative formulations '.... Lead in principle to nonsymmetric 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 hyperelastic stressstrain 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 strainspace 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 elastoplasticity based on the multiplicative decomposition. Despite their success, these approaches involve modifications, and often a complete reformulation, of the standard closestpoint algorithms of the small strain theory. From a practical standpoint the implication is that the implementation of classical models needs to be considered on a casebycase basis in the large strain regime. In a later study, Simo [20] proposed a stateoftheart 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 closestpoint projection algorithm of any classical simplesurface or multisurface 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 closestpoint 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 timestepping 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 abovementioned theory [20] of multiplicative decomposition to derive explicit expressions of consistent tangent moduli for the proposed elastoplastic 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 (nphase) 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 noninteracting and nonpolar. 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 singlephase 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 twophase 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 volumeaverage 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 PaVznzna PaVxAxn 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 PaVxAxn 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 [2126] 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 nonpolar 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) afn lP va.a d + r (lpava.v va .n)dr, Dat CaHa df 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 qIrqa 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 positivevalued 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 nphase continuum considering n = 2 (solid and fluid phases). Following the derivations in Section 3.4, necessary balance laws for soilwater 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 soilwater mixture as p+div(p)= 0, (3.40) at where i is the volumeaverage 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 fcndrf= (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 (nonsymmetric) first PiolaKirchhoff 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. Ft = p +Pw, (3.43) where Ps = Joa Ft and P" = Joa Ft are the first PiolaKirchhoffpartial 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 PiolaKirchhoff 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 pullback mass density of the soilwater 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 massdensity 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 soilwater 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 nonmechanical 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 lefthand 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 righthand side of (3.49) as J(pag+aa)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 soilwater 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 PiolaKirchhoff total stress tensor such as S= F P = F Ft = JF1' iFt. (3.57) Following (3.56) and (3.57), additive decomposition of S can be written as S= S+9C1, (3.58) where S is the second PiolaKirchhoff 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 Jn )tw =rw= :dw SW . 0 2 Balance of energy for the soilwater 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 soilwater 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 soilwater 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 nonmechanical 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 soilwater 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 soilwater 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 Jnw )lw JY = (3.67) Jp From (3.62), (3.65), and (3.67) localized version of entropy inequality for soilwater 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 welldefined 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 wellposed 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. Pushforward 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 pullback 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: PpoTig)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 ipoT. 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 semidiscretization 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 kstep 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+im + Ati3mfn+lm)= 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 [5861]. These kstep 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+lm)+ AtPofn+l = 0 (4.22) m=0 with k + 1 unknown coefficients which one can choose by forcing a kstep 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+lm) Atlofn+l =0. (4.23) m=1 , See [5] for determination of coefficients for BDFk scheme. (4.23) can be made more generalized by introducing onestep recurrence relation for f such that dn+i ((amdn+im) AtPo[lfn+i+(1P)fn = 0. (4.24) O < p 1 is the parameter of generalized trapezoidal method. Some of the wellknown families of 1methods are presented in the following. Table 4.1 Generalized trapezoidal method P Method 0 Forward Euler 1/2 CrankNicolson 1 Backward Difference Unconditional stability is achieved for any At if 1 > 1/2. In general, the nonlinear responses of interest are dominated by lowfrequency 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 highfrequency components but at the expense of accuracy. On the other hand, CrankNicholson 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 timeintegrated variational forms of HAt (, I, w) given as HAt Hn.W) =v/ f Jn+i SamJn+im dV P0 oJb(grad JV)n+ +(1 )(gradl .J)n] dV (4.25) Po10BPQn+l +(lP)Qn]dA= 0, HAt(,Hw)= Jn+l mamJn+m dV PoJg0 p(GRADyV)++(1P)(GRADw.)n]dV (4.26) Poj oIPQn+ +(1P)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 platelike 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 nonmechanical 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 Cdr 0 (4.28) dt 2 dt where d = Sym(l) is the rate of deformation tensor, I = F FI 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 CauchyGreen tensor, b', provided that b* satisfies an objective transformation. An elastoplastic 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 (stresslike) internal variables characterizing the hardening response of the soil. t e R" is a vector of internal plastic variable (strainlike) 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 timederivative 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 CauchyGreen 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 elastoplastic 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 secondorder 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+LBgradk .(Pgrad9 +JgdVJ lQdA. (4.53) P Pwg CHAPTER 5 HYPERELASTICPLASTIC MCC MODEL 5.1 Introduction Elastoplastic 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 CamClay (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 elastoplastic 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 elastoplastic 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 nonconservative 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 CamClay 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 twoinvariant stored energy functions [68] is employed in this study which includes pressuredependent 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 [6971]. 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 ==P8; 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 timevariation 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. Twoinvariant yield function of MCC model [28] is given by the ellipsoid 2 f= f(p,q, pc)= +P(PPc)=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 pq plane. Yield Surface Figure 5.1 Yield surface of the MCC Model in pq 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 elastoplastic 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= ls'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) + rI2q  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 straingradient 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 =[dlb22cd2b21c+det(b)Ktr , e O p (5.49) a2 = 1 (d2b Idib12 , 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 2de2t(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 elastoplastic tangent moduli aep =D ep s s + 2q (Ii@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 elastoplastic (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 elastoplastic 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 CamClay model. CHAPTER 6 HYPERELASTICVISCOPLASTIC MCC MODEL 6.1 Introduction Hyperelasticplastic MCC model is further extended for viscoplasticity to model timedependent 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 elastoviscoplastic constitutive models have been proposed to describe the theological behavior of clay. Most elastoviscoplastic constitutive models can be classified as overstress models or nonstationary flow surface models. Overstress elastoviscoplastic 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], BaladiRohani 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 nonstationary 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 nonstationary 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 timerate flow rule offers the capability of simulating timedependent 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 Aaf; 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 = 2II3e 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 elastoviscoplastic 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) + 3tr2q (I5@8n . 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 elastoviscoplastic tangent moduli as follows: aeP = @ii+ iDep@8 Ofi+ J3 12 (6.21) + 4(In@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"' F1 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 soilwater 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 +OFt)+ Jpg. (7.6) Taking the variation yields E = DIV(SP + Ft + 0 F)+ 5(Jp)g BE = DI (7.7) =DIV(A :F +0Ft +SOFt)+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(F1 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(JF1 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(OFt GRADtSu Ft) +DIV(68Ft)+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 PiolaKirchhoff effective stress tensor S as 3 S=F' Ft = XPAM(A); M(A) =F1.m(A) Ft. (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 pushforwards 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 pushforward 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 CauchyGreen tensor b having distinct eigenvalues )X, 2 and 1 d(A) takes the form d(A) = [Ib b@b+I3A2(l1I) 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 = JF1 v)= 5JF1 v + JF1 v+JF1 *v. (7.22) From (7.4a), (7.9) and (7.3b), one can express 8JFiv=DIV(JF1 8u)1 (7.23) ( '(7.23) JF1 v =JF1 grad u v = JF1.GRAD8u F1 v. Substituting (7.23), one can rewrite (7.22) as V = 5(F1 v)= DIVJFI1 .u vJF1 GRADBu F + JF1 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 = F1 k Ft is the pullback 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(JF1 .uFt + JGRADt8u. (7.28) Substituting (7.28) and rearranging, (7.27) can be rewritten as 8V = K GRAD5 + [DIJF1 8u)Ft + JGRADt5u 1 (7.29) SK (GRAD0 +J.IF G 8Pw t where 6K=F .{2Sym(GRADu KFt)(I+eo)DIV(JF' .u)ak}.Ft.(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 soilwater 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 F1). + 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 (8div08gradti: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 straightforward 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: PpoT g)dV 5(Ti.t)dA =G+J 8(GRAD q:P)dV+J 8(GRAD : (OFt))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: (bFt))= GRADI :Ft 8 GRAD : (Ft GRADtSu Ft) Upon substitution of (7.10), variation of the third integrand yields 8(poT g)= PwDIV(JF1I 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:FtdV IB RAD: t(.GRADt8uFt)dV (7.40) f pwDIV(JF1 8u gdV T8tdA. (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 UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EZZBLI29T_5JFZN0 INGEST_TIME 20130408T23:19:33Z PACKAGE AA00013546_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES xml version 1.0 encoding UTF8 REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID ESRBCC5EZ_GZSBXH INGEST_TIME 20130312T14:21:24Z PACKAGE AA00013546_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 