2 Thermodynamics-based constitutive framework for unsaturated soils 2

Li,X.S.(2007).Ge´otechnique 57,No.5,423–435423Thermodynamics-based constitutive framework for unsaturated soils.2:A basic triaxial modelX.S.LI*Following the theoretical framework established in Part 1,a basic elastoplastic unsaturated soil model is formu-lated.The model is a combination of a soil skeleton model and an SWCC (soil-water characteristic curve)model.For simplicity,the model deals only with isotropic soils in triaxial compression space under a constant pore air pressure.The paper shows that the central element in the model that differentiates itself from a saturated model is its dual plastic hardening,which not only reflects the evolution of the internal structure of the soil skeleton but is also associated with the evolution of the solid–fluid interactive structure.For more flexible shear stress–strain responses,the non-associated flow rule (the alpha model)by Collins &Kelly,based on a thermody-namic analysis,is adopted.The present model introduces six new material constants,which can be calibrated following a systematic procedure.The model successfully simulates the multiple series laboratory tests reported by Sivakumar &Wheeler,including wetting,consolidation and shearing,on a compacted clay.KEYWORDS:constitutive relations;partial saturation;plasti-city;pore pressures;stress path;suctionUn mode`le e ´lastoplastique de sol non sature ´est formule ´sur la base du cadre the´orique de ´fini dans la partie 1.Il combine un mode`le de squelette solide et un mode `le de courbe caracte´ristique sol-eau (SWCC).Pour simplifier,ce mode`le ne traite que de sols isotropes dans un espace de compression triaxiale sous une pression d’air dans lespores constante.Cet article montre que l’e´le ´ment central du mode`le qui se distingue d’un mode `le sature ´est son double durcissement plastique,qui refle`te non seulement l’e ´volution de la structure interne du squelette solidemais correspond e´galement a `celle de la structure inter-active solide-fluide.On adopte la loi d’e´coulement non associe´e de Collins et Kelly (le mode `le alpha),base ´e sur une analyse thermodynamique,pour obtenir des re´ponses plus flexibles de contrainte-de´formation au cisaillement.Le mode`le pre ´sente ´ici introduit six nouvelles constantes de mate´riau,qui peuvent e ˆtre e ´talonne ´es suivant une proce´dure syste ´matique.Il simule avec succe `s les essais de se ´ries multiples rapporte ´s par Sivakumar et Wheeler,re´alise ´s en laboratoire sur une argile compacte ´e,incluant la mouillabilite´,la consolidation et le cisaillement.INTRODUCTIONA thermodynamics-based unsaturated soil modelling frame-work has been presented in the preceding companion paper (Li,2007),referred to as Part 1hereafter.The objective of the present paper (Part 2)is to demonstrate a specific model constructed in accordance with that framework.For effectively presenting the main idea,only a basic model is described in the sense that it deals only with isotropic soils in triaxial compression space with a constant pore air pressure.Extending the model to general 3D stress–strain space with a variable pore air pressure is straightforward;and modelling anisotropic features (includ-ing cyclic shearing)is a general,important,but separate issue beyond the scope of the present work,and therefore will not be addressed in this paper.In this model,because the pore air pressure is assumed to be constant,the behaviour of an unsaturated soil is described by two interrelated models:one for the solid skeleton,and the other for the pore fluids (referring to equations (88a)and (88b)in Part 1).As shown in Part 1,the former is constituted through specifically defined elastic moduli,plas-tic yielding,flow and hardening,and the latter is a model of soil-water characteristic curves (SWCC).All these model ingredients and their integration will be addressed in depth,supplemented with necessary references.The calibrationprocedure of model constants will also be presented in detail.To examine the performance of the model,the test results reported by Sivakumar &Wheeler (2000)and Wheeler &Sivakumar (2000)on compacted speswhite kaolin were simulated.A comparison between the simulated results and the test data shows that the model broadly captures the essential behaviour of the tested soil through all the stages,including wetting,consolidation and shearing,of the tests,indicating that the modelling framework proposed in Part 1and the specific model presented in Part 2(this paper)indeed have a sound basis.MODEL OF SOIL SKELETONIn triaxial compression space,the stress and strain tensors are reduced to two-component vectors~pq &'¼~óa þ2~ór ðÞ=3~óa À~ór &'(1a)andåv åq&'¼åa þ2år 2åa Àår ðÞ=3&'(1b)where åv and åq are the volumetric and deviatoric strainsrespectively,and ~pand q are the mean normal and deviato-ric stresses conjugate to d åv and d åq respectively in work done on the soil skeleton.In the definition,åa and år are theaxial and lateral strains respectively,and ~óa and ~ór are their corresponding stress conjugates.Manuscript received 11August 2005;revised manuscript accepted 25January 2007.Discussion on this paper closes on 3December 2007,for further details see p.ii.*Department of Civil Engineering,Hong Kong University of Science and Technology.Elastic relationsThe elastic response of the soil skeleton is described by the commonly used incremental relationsd~p¼K dåev ¼v0~pkdåev(2a)d q¼3G dåe q¼33K1À2íðÞ"#dåe q(2b) where K and G are the elastic bulk and shear moduli respectively;í(Poisson’s ratio)and k(a Cam clay para-meter)are two material constants;and v0is a reference specific volume at which the volumetric strainåv¼0.Note that the shear modulus so defined depends on the mean normal stress~p,making equation(2b)hypoelastic:this means that the relationship given by equation(2b)is not exactly elastic,in the sense that recovery of the devia-toric stress upon an arbitrary closed cycle of strainåe q is not guaranteed.The impact of this effect on the overall stress–strain response,however,is considered negligible,because of the dominance of plastic deformation.In other words,in this model the elastic and plastic deviatoric strain increments aredenoted as dåeq and dåpqinstead of general infinitesimalsäåeqandäåp q(referring to the notation in Part1)only in an approximate sense.On the other hand,equation(2a)is an integrable equation,which makes the relationship between~pandåev exactly elastic(path independent).Although this factis widely known from study of the Cam clay models,it is worth reiterating that,because of this exactly elastic relation-ship,the effect of the plastic volumetric strainåp v is guaran-teed to be path independent,and consequently it can be treated as a perfect internal state variable in formulating the hardening laws.In the present model,the internal structure of the soil skeleton,H9s,is uniquely characterised byåp v. Yield function and plasticflowIn the triaxial space,the dissipative equilibrium equation (equation(37)in Part1)is reduced top D Dþq DÀv s¼~pÀpÆðÞDþqÀqÆðÞÀr s iæs i¼0(3)where p D and q D are the dissipative stresses;D¼dåp v=dåp q is the dilatancy such that the direction of the plastic strain increment m is equal to f D1g T and the loading index L s is equal to dåp q;and pÆand qÆare the back mean normal and deviatoric stresses.In this model,as the material is treated as isotropic,the yield function must be unbiased with respect to all three principal stress axes,and accordingly the deviatoric back stress qÆmust vanish identically,resulting in q D q.By setting q D¼q,takingåp v andåp q as two internal variables,and then following the idea of Collins&Kelly (2002)concerning a non-associatedflow rule(the alpha model),the above equation is converted into the more concrete formp D Dþq¼ð0DffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiD2þM2R2p!|fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}æs1D|{z}r s1þð0M2R2ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiD2þM2R2p"#|fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}æs21|{z}r s2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}v s(4)where M is the well-known critical state stress ratio,andð0 and R are state functions.Equation(4)indicatesp D¼ð0DffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiD2þM2R2p(5a) andqMR¼ð0MRffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiD2þM2R2p(5b)Adding the squares of the above two equations together results in the yield function in the dissipative stress space,f s¼M2R2p2Dþq2ÀM2R2ð20¼0(6)Equation(6)represents an ellipse in dissipative stress space. By denoting the values of~p on the yield surface intercepting the hydrostatic axis(q¼0)as p b and p0+p b respectively, one hasð0¼p0=2¼pÆÀp b(7) The dilatancy is then obtained by imposing the normality rule in the dissipative stress space:D¼@f s=@p D@f s=@q¼M2R2p Dq¼M2R2~pÀpÆðÞq(8)where R is given,based on the original proposition of Collins&Kelly(2002),asR¼Æþ1ÀÆðÞp Dþð0ð0¼Æþ1ÀÆðÞ~pÀpÆþð0ð0(9) in whichÆis a material constant.In the present model,ð0¼ð0(åp v,ns,S r)and pƼpÆ(åp v, ns,S r).Such functional dependences make the thermody-namic forcesæs1andæs2depend not only on the internal state variable of the soil skeletonåp v but also on thefluid state variables ns and S r(referring to equations(16)and(17)in Part1).Combining equations(6),(7)and(9)transforms the yield function into the true stress space,reading asf s¼f s~p,q,p0åp v,ns,S rÀÁ,p båp v,ns,S rÀÁÂüf s€p,q,p0ðÞ¼M241ÀÆðÞ2€p4p2À41À3Æþ2Æ2ðÞ€p3p0ÀÆ4À5ÆðÞ€p2ÀÆ2€p p0!þq2¼0(10)where€p¼~pÀp b(11) is the shifted mean normal stress.Fig.1shows the yield function plotted in true~pÀq stress space.It can be seen that p b is the back mean normal stress(referring to equation (78)in Part1)that controls the position of the yield surface in the true stress space,and p0characterises the size of the yield surface.Note that becauseð0¼ð0(åp v,ns,S r),p b and p0are,in general,both functions of the state variablesåp v,ns and S r.The condition of consistency of the yield function reads asd f s¼@f s@~pd~pþ@f s@qd qþ@f s@p0d p0þ@f s@p bd p b¼@f s@€pd€pþ@f s@qd qþ@f s@p0d p0¼0(12)424LIwith the partial derivatives given by @f s @~p ¼À@f s @p b ¼@f s@€p¼M2"161ÀÆðÞ2€p 3p 20À121À3Æþ2Æ2ðÞ€p 2p 0À2Æ4À5ÆðÞ€p ÀÆ2p 0#(13a)@f s @q ¼@f s@q ¼2q (13b)@f s @p 0¼@f s@p 0¼M 2À81ÀÆðÞ2€p 4p 30þ41À3Æþ2Æ2ðÞ€p3p 20ÀÆ2€p "#(13c)It can be seen that ì¼{@f s /@~p @f s /@q }T is in general notco-directional with m ¼{D 1}T ,indicating that a non-associative plastic flow is materialised.As demonstrated by Collins &Kelly (2002),the above yield function (with p b ¼0)is reduced to the modified Cam clay yield function when the parameter Ƽ1,and only in this special case is the associated flow rule resumed.It should be remarked that the yield function,as well as the corresponding flow rule presented above,was originally proposed by Collins &Kelly (2002)based on their thermo-dynamic framework for saturated soils.The adoption of this yield function is based on the consideration that its non-associated flow rule offers more flexibility in modelling shear stress–strain behaviour by adding only one more material constant,as compared with the modified Cam clay yield function.However,the latter,as well as many other existing yield functions,can also be adopted.Hardening and plastic moduliThe present model includes both isotropic hardening and kinematic hardening,characterised by the internal state vari-ables p 0and p b respectively.It is noted that in both the Cam clay and the Collins &Kelly models,p b 0,and p 0is assumed to depend only on the plastic volumetric strain.This is valid only for saturated soils,as it involves no impact of unsaturation.In the present unsaturated soil mod-el,as attention was called to earlier,the two hardeningvariables p b and p 0depend not only on åp v but also on ns and S r ,the state variables of the pore fluids.Such depen-dences are consistent with the theoretical framework pre-sented in the companion paper (Part 1),and are vital for a proper modelling of essential unsaturated soil behaviour.The hardening response is investigated here through iso-tropic compression.Fig.2shows a simple mechanism of aparticle sliding under macroscopically isotropic stress ~p.It is assumed that an irreversible change in volume is attrib-uted mainly to particle sliding (although the macroscopic shear strain may be zero),which repacks the particle assem-blage into a different density state.During isotropic com-pression,both the normal and shear forces at particlecontacts are associated with ~p.When ~p Àp b is less than a threshold p 0,the contact shear force is,on average,less than the product of the contact normal force and the coefficient of friction between particle minerals,such that sliding is notmaterialised.When ~pÀp b ¼p 0,the resultant normal force times the coefficient of friction Øis equal to the resultant shear force,as illustrated in Fig.2.This is the macroscopic yield condition.By scaling the resultant shear force to Øp 0,the resulting normal force is equal to p 0,which is a state function that can be written in the formp 0åp v ,ns ,S r ÀÁ¼p ååp v ÀÁþp ìåp v ,ns ,S r ÀÁ(14)where p åis a parameter of the soil skeleton solely related to its internal structure (H 9s ),and p ìreflects the impact of the solid–fluid interaction (H 0s ).Note that for a soil in saturated or dry states H 0s disappears,and so does p ì(to be consistent with saturated models,in saturated or dry states the back-pressure p b ¼0as well).This indicates that p åis simply the yield confining pressure for the soil in saturated or dry states,given by the well-known empirical observation (Schofield &Wroth,1968)p ååp v ÀÁ¼eN Àv p ºÀk ¼e N Àv p 0þv 0åp v ºÀk (15)where º,k and N are the well-known Cam clay constants,denoting the slopes of the normal consolidation line and therebound lines in the v Àln ~pplane and the value of the specific volume on the isotropic consolidation line at ~p¼1(note that in saturated states p b vanishes,so it plays no role)respectively;and v p and v p 0are the plastic parts of the specific volume and its initial value respectively.As shown in the above equation,p åis uniquely related to v p ,meaning that p åis an internal state variable of the soil skeleton.That is,p åevolves only when the soil skeleton experiences an internal structural change (an irrecoverable volume change).In this basic model,it is further assumed that p å(or,equivalently,v p or åp v )is the only internal state variable of the soilskeleton.1020304050p : kPa: k P aq ~Fig.1.Yield surface in ~p2qplane p ns S p p ns Sp p p ns S p ns S 0p v r p v pv r p vp vrbp vr(,,)()(,,)()(,,)(,,)εεεεεεϭϩϭϩϩεµεµ~Fig.2.Sliding block model for yielding under isotropic stressconditionTHERMODYNAMICS-BASED CONSTITUTIVE FRAMEWORK FOR UNSATURATED SOILS 2425On the other hand,p ìcharacterises the impact of solid–fluid interaction (meniscus effects)on the soil skeleton,which is related not only to åp v but also to the fluid state variables ns and S r .This means that p ì,and consequently the size of the yield surface (p 0),may change even if the internal structure of the soil skeleton remains unchanged.This is the case where äH 9s ¼0but äH s ¼0(referring to equation (43)in Part 1).For example,during wetting,as suction becomes smaller,p ìis reduced,as is p 0.However,itis entirely possible that ~pÀp b is smaller than p 0such that the soil skeleton still remains in elastic states until ~pÀp b ¼p 0.As p ìcharacterises the impact of meniscus effects on the soil skeleton,it is zero when ns ¼0(soil is saturated)or S r ¼0(soil is dry).For convenience,an intermediate para-meter ì(åp v ,ns ,S r )>1is introduced,with ì¼1for ns ¼0or S r ¼0,such that p ì¼1Àr ðÞìc À1ðÞp å(16a)orp 0¼p ìþp å¼ìc Àr ìc À1ðÞ½ p å(16b)where r and c are material constants.In this model,ìis given by the following empirical relation (other specific functions may be proposed when more test data become available):ì¼ìåp v ,ns ,S r ÀÁ¼1þap 1Àb atm nsS r ðÞb p å¼1þap 1Àbatm nsS r ðÞb ev pÀv 0åp v ÀN ºÀk(17)where a and b are two additional material constants,and p atm is the atmospheric pressure,serving as a reference to make the material constants dimensionless.The parameters a ,b ,c and r are all the material constants added to conven-tional saturated model constants to characterise the behav-iour of the soil skeleton.The calibration procedure of these new constants will be presented later.It can be seen that equation (17)satisfies the constraint that ì¼1for ns ¼0or S r ¼0;and ìtends toward unity as soil becomes denser,reflecting the general trend that soil becomes more stable (less affected by solid–fluid interaction)when it has a denser structure (Vicol,1990).As p ìand the back mean normal stress p b both originate from meniscus effects,they are homogeneously related.As a first approximation,p b is assumed to be proportional to p ì,asp b ¼r1Àr p ì¼r ìc À1ðÞp å(18)such that the isotropic yield stress is given by the simple relation~p¼p 0þp b ¼p åìc (19)It is interesting to note that both ~p¼þsS r and p 0+p b are functions of sS r ,with ~p<p 0þp b always.If initially ~p,p 0þp b (say,immediately after compaction),when wet-ting proceeds under a constant net pressure p ,sS r decreases,and so do ~pand p 0+p b .If the decreasing rate of p 0+p b is greater than that of ~p,at a particular value of sS r ,~p may be equal to p 0+p b ,signalling the onset of plastic yielding.Any further decrease of sS r after that instant will be associated with a permanent volume reduction.This is the case of collapse.As illustrated in Fig.3,collapse may occur only when p is greater than the initial p å,indicating that,during wetting,collapse is likely to occur if the applied pressure (p )is high and/or the material is loose (p åis low).The evolution of the isotropic hardening state variable p 0is given by differentiating equation (16b),with ìreplacedby equation (17),as follows:d p 0¼1Àr ðÞìc þr ÂÃd p åþc 1Àr ðÞp åìc À1d ì¼1Àr ðÞìc þr ÂÃ@p å@åp v þc 1Àr ðÞp åìc À1@ì@åp v &'d åp vþc 1Àr ðÞp åìc À1@ì@S r þ@ì@ns ðÞ@ns ðÞ@S r!d S r ¼v 0p 0þp b ðÞºÀk p 00b À1Àr ðÞc ìÀ1ðÞì!d åp vþ1Àr ðÞbc p 0þp b ðÞìÀ1ðÞns þS r ÇðÞnsS r ìd S r(20)and the evolution of the kinematic hardening state variable p b is given by differentiating equation (18),with ìreplaced by equation (17),asd p b ¼r ìc À1ðÞd p åþcr ìc À1p åd ì¼v 0p 0þp b ðÞºÀkp b p 0þp b Àr c ìÀ1ðÞì !d åp vþrbc p 0þp b ðÞìÀ1ðÞns þS r ÇðÞnsS r ìd S r(21)where the relations given in equations (15)and (17)have been invoked.In the above equation,Ǽ@(ns )=@S r is to be determined by a separate SWCC model (referring to equa-tion (69)in Part 1),which in the present work features the hysteretic response of the scanning curves (see next section)so that Çis hydraulic path dependent.By combining equations (12),(20)and (21),invoking thedefinition d åp v ¼D d åpq ¼D h L s i ,and then referring to equa-tion (45)in Part 1,one obtains the plastic moduliK 9p ¼v 0p 0þp b ðÞD ºÀk (@f s @€p p b p 0þp b Àrc ìÀ1ðÞì!À@f s @p 0p 0p 0þp b À1Àr ðÞc ìÀ1ðÞì !)(22)andInitial p sS rFig.3.Illustration of evolution of p 0+p b during wetting426LIK 0p¼@f s @€p r À@f s@p 01Àr ðÞ"#bc p 0þp b ðÞìÀ1ðÞns þS r ÇðÞnsS r !(23)where the partial derivatives @f s =@€p and @f s =@p 0havebeen given in equations (13a)and (13c),and the dilatancy Dhas been given by equation (8),in combination with equa-tions (7)and (9).Note that when yielding takes place under isotropic stress conditions,€p ¼p 0so that @f s =@€p ¼À@f s =@p 0(refer to equations (13a)and 13(b));and,accord-ingly,the roles of r in equations (22)and (23)are elimi-nated,which means that the material constant r plays no part in isotropic consolidation and/or isotropic wetting/dry-ing.Note also that as Çis wetting/drying path dependent,K 0p is wetting/drying path dependent,too,which means that the soil deformation induced by hydraulic actions depends on the history of the hydraulic process.MODEL OF HYDRAULIC RESPONSEThe bounding surface SWCC model proposed by Li (2005)based on phenomenological observations on s –S r scanning curves is used in the present study with a minor modification:that is,the soil-water characteristic curves are expressed as relations between ns and S r ,rather than be-tween s and S r ,because ns and S r are conjugate in volume-specific work,as shown in equation (19b),Part 1.The model assumes two primary soil-water characteristic curves for a given soil.One is f d (ns ,S r )¼0for primary drying (presumably starting from a full saturation state)and the other is f w (ns ,S r )¼0for primary wetting (presumably starting from a completely dry state),where s is called the image suction,referred to as the values of suction on the primary curves.These two curves serve as the bounds of all scanning curves.In addition,a weighted back suction (referring to equation (59)in Part 1)is also defined,which is the product of porosity n and the value of suction at which the last reversal of the hydraulic process from wetting to drying,or vice versa,took place.Because the weighted back suction changes discontinuously (only when the hy-draulic process changes the direction),the value of is piecewise (within each scanning curve)constant.In the model,the wetting/drying state is specified by the degree of saturation S r and three weighted suction values,ns ,and ns ,with representing the impact of the wetting/drying history.The relations among these state variables and their increments are illustrated in Fig.4.Refer to Li (2005)for more details (with slightly different notations:and d ¼Æand â,respectively,in Li,2005).The general incremental relationship among ns ,S r and åhas been given by equation (69)in Part 1.In the present model,partially because of lack of systematic experimental results,the relation is simplified by assuming ˛¼0.That is,the impact of solid–fluid interaction on SWCC is consid-ered to be negligible.Therefore d ns ðÞ¼H 9p M f H 9p þsgn äS p rÀÁ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}Çd S r ¼Çd S r (24)Following Li (2005),the stiffness Çis given byǼj ln ns =âÀÁj dH j ln ns =âÀÁj d þj ln ns =âÀÁj d H À1ðÞ"#s s d ns ðÞd S r (25)where H and d are two model parameters,which can be either constants or state functions.In the present work they are both treated as constants for simplicity.The parameter d controls the ‘fatness’of the hysteresis loops formed by the scanning curves.It can be seen from equations (24)and (25)that,when a scanning curve starts right after a change in wetting/drying direction,ns ¼â,and the ratio d ln (ns )=d ln (ns )¼H ,which controls the initial slope of the scan-ning curve;and as ns moves on its way towards ns ,the slope of the curve varies smoothly,and when ns gets close to ns ,d ln (ns )=d S r %d ln (ns )=d S r ,indicating that a scan-ning curve is an asymptote of its respective primary SWCC.The model described above can be used in conjunction with any propositions for the primary wetting and drying curves.In this study,a version of the SWCC equation proposed by Fredlund &Xing (1994)is used,which is S r ¼1ln exp (1)þns kp atml "#()m(26)where k ,l and m are curve-fitting parameters.In this model,there are two sets of such constants,k d ,l d ,m d and k w ,l w ,m w for the primary drying and wetting curves respectively.STIFFNESS MATRIXAs stated earlier,for simplicity,the pore air pressure in the present model has been assumed as a constant (pore air spaces connect to a constant pressure source,say,the atmo-spheric pressure;and loading is sufficiently slow),so that equation (88)in Part 1consists of only the first two simul-taneous equations,which in the triaxial space can be written asd ~pd q d ns ðÞ8<:9=;¼¸11¸12—1¸21¸22—20Ç2435d åv d åqd S r8<:9=;(27)where ˛has already been set at 0.In the above equation,the coefficients ¸ij and —i are the components of the stiffness matrices ¸and —given by equations (53)and (54)of Part 1,which in the triaxial space are writtenas0·20·40·60·81·0ns : kPaS rFig.4.Illustration of bounding surface SWCC modelTHERMODYNAMICS-BASED CONSTITUTIVE FRAMEWORK FOR UNSATURATED SOILS 2427。