Acta Polytechnica Vol. 43 No. l/2OOg Numerical Impleme ntation of rsotropic Consolidation of Clayey Soils T.Janda, P. Kuklft, M. Sejnoha 1!x paler llforts on implernentation of seaeral numerical _techniques to solae a set _of goueming equations resutting frorn simple onedimensionalisotropit mnsol'i'dation of soik thatbehaae according to.i\ cay.ck2 mode.l."Tiree dffiTenimahods of ntolif tlre rq:*r:i"", aconsol'idat'ion, namcly the collocation ry,et!,od, th1 f;nite ao.lume lnethod and the finite element neiiod,, orc piu*nte6. Apart from eaaluatrng their effrciency, whlch becornes particul'a,rly crucinl uhen-implementing these teihniques in the framework of an optimization problem aimedaj hmi!'g thc model parameter\ 1' set of parameters of a Cam clal"modcl d,riaing the time'dzpenfunt ,"tponiu of the soiis (deformation dcpendcnt aariation of the cofficiznt of permeability and, preconsoiidation pressure) is ako d,iscussed,. Keywords: Cam clay, coeffi'cient of permeability, consoli.dation, finite aolume method, aoid. ratio, genetic algorithm. I Introduction Owing to its relative simpliciry isotropic consolidation is often viewed as a basic tool for inferring the material parame- ters of various constitutive models. While the respective con- stitutive model is usually well defined and tested for the de- scription of a certain type of soil, the values of the associated material parameters are mostly lacking, leading to an exten- sive laboratory program to determine them, On certain occa- sions, howeve4 this problem can be confined to a single labo- ratory test combined with numerical simulations. In particu- lal this paper advocates the use of a simple isotropic consoli- dation test to infer the basic material parameters describing critical state models such as the modifred Cam clay model. Among others, the modified Cam clay model has been often the choice for constitutive models for a realistic repre- sentation of the inelastic behavior of clayey soils, particularly when the deformation of the solid phase is of the main concern. Its selection has been promoted by the ability of the model (with only minor modifications) to capture a number of important phenomena associated with the soil-water inter- action. Wth reference to consolidation, the need was long ago highlighted for at least a bilinear form of the consolidation line to account for the prior loading history (unloading from a certain level of preconsolidation pressure prior to subse- quent loading), which signif,rcantly affects the shape of the pressure dissipation cuwe. Successful implementation of the model further requires a deformation dependent formulation of the coefficient of permeabiliry. When expressed as a function of the actual void ratio this can significantly improve the material response, particularly when the degree of consolidation increases [8]. See also [4, 9] and Section 4 for further discussion. The above requirements then give rise to five basic material parameters to feed into the constitutive model. These parameters are summarized in Section 2, which briefly reviews the theoretical formulation of isotropic consolidation. Manual tuning of these parameters is usually not very effi- cient, and should be combined with a suitable and reliable optimization technique to determine their values [0]. Here the efficiency of the selected numerical technique for solving the governing equations is of paramount importance. There- fore, a critical evaluation of the efliciency of these methods, discussed in Section 3, is one of the goals of this paper. 8 2 Governing equations Referring to the experimental measurements carried out in the triaxial apparatus, isotropic consolidation can be view- ed as a one phase ow in a fully saturated deforming medium undergoing small deformations. Neglecting the body forces, the hydrostatic state of srress maintained during the experi- mental measurement gives or(x,y,z,t) = o r(x,y,z,) = o,(x, y,z,t) = c*(t), ( I ) where o, is the total mean stress. Following the Terzaghi- -Fillunger concept of effective srresses this quantity can be expressed in terms of thepore pressure p' and the effective stresses benveen grains offI as -cFo*=67;'-p'. (2) Assuming full saturation (S,= l, [6]) the pore pressurep' equals the pressure in the liquid phase p'. Referring to the experimental conditions, the total mean stress remains con- stant throughout the consolidation process. The assumed stress homogeneiry together with Eq. (2) then provides o*=of -i'=o (3) where ( ' ) represents the time derivative a( ) I at . Thansport of the liquid phase throughout the soil sample can be described by the following set of equations: Transport eEtation r =-E{gradpu, (4)J ^,wI where ,fl is the mass flux of pore water, y' = gp' is the spe- cific weight of water, p' is the intrinsic mass density, and K represents an instantaneous coefficient of permeability. Balance equnfion reads P- i, t divrl? =9. (5) The volumetric strain e.. follows from the Constitutiae eEntiun e(t)-Vo K ,€r,=i1-a=-; :'_ h(-o;tr(rD, ":f , ":f; (6) l+eg I+eg e(t) - eo l"e,=ia-l=-;:tn(-of;1r), of .eff Q) I-t €6 l-f 01y Acta Polytechnica Vol. 43 No. l/2003 (13) {-6'# ) h(-of ) Fig. I : Bilinear form of the consolidation line derived in the case of the modified Cam clay model llom the bilinear consolidation line, Fig. l. The initial branch, often referred to as the r - line, gives evidence of the previous stress history and represents the effect ofoverconsolidation. The slope discontinuity between the k and l, - lines can be identified with the structural strength of soil given in terms of a certain level of the effective mean stress ofitr (pre- consolidation pressure). Differentiating Eq. (7) with respect to time gives the rate of volumetric strain in the form ., /.\ i'(t) 7, b'U) crtt", - r+ es * ^;nt' (8) Substituting Eqs. (a) and (8) into Eq. (5) and taking into account the actual triaxial set-up, in which only the bottom face of the cylinder is drained, leads to -+";ut,l 9 lrtry ael(rl l - b' (t) = o"l*)\ o.t oz ) It has been verified experimentally that in the case of iso- tropic consolidation a simple power law written as [8] K(t\ (e(i\'',)', =l -.-, | , (10)K6 \to ) represents the soil behavior fairly well. The dependence of the actual void ratio r on the e{fective mean stress, Eq. (7), to- gether with Eq. (10) provides aKQ) = mK(A a?'(t). (')az eQ)c'f Q) az Introducing Eq. (11) into Eq. (9) finally yields b.(t\ =- K(r)(l + 4o) . vt)t 6(L 5 E @o E 0) I 00 I 1 (e) lnfryl ."*u'fff (r2) A similar equation can be derived for the unloading branch when replacing l. by rc and ee by as in Eq. (12) In(-off ) lreWl ."*u'affi) 3 Numerical implementation The purpose of this section is to evaluate several numeri- cal methods commonly used when solving the consolidation problem. As mentioned above, the main objective remains computational efliciency of individual numerical techniques. Frrst, we recall the collocation method successfully imple- mented in [8]. The implicit finite control volume scheme combined with the finite difference method is explored next [6]. Frnally, our attention is given to the finite element method, which has proved to be an efficient tool particularly when taking the volume changes of the porous skeleton into account [2, 9]. Before proceeding with the present formulation of the individual numerical techniques, we recall awell known draw- back of standard formulation attributed to the point of slope discontinuity along the bilinear consolidation line. This phe- nomenon, evident fiom Fig. 2 (see also [7]), results in a typical eo = 0.5, Ko = 2.0'10-to[ms-l r = 0.01, l, = 0.'l -0 1000 2000 3000 4000 Time [min] Fig. 2: Time variation of pore pressure with no smoothing of the bilinear consolidation line double S type of curve where individual segments of this curve are linked to a given slope of the bilinear consolidation line. This, howeve4 is in contradictionwith actual experimen- tal measurements. To avoid numerical difliculties around the point of discon- tinuity we equipped the original diagram in Fig. I with a cubic parabola to smooth out the transition zone. The numerical integration ofEqs. (12) and (13) is then carried out under cer- tain simplifying assumptions. In particula6 we first introduce an instantaneous modulus l" with corresponding initial void ratio 0g displayed in Fig. 3. These parameters, evaluated at the middle of the current time step, are then introduced in Eqs. (12) and (13), thus replacing the original variables r, 1", eo and eg, respectively. Furthermore, the pore pressure de- pendent material parameters K and e are derived fiom an ini- tial prediction ofpore pressure taken at the end ofthe current time step. Details are given in Sections 3.1-3.3. Acta Polytechnica Vol. 43 No. l/2005 b' (t 1 a) = -u (" ,D-S :!-*l,r* o)' tn1-oif ) In(-o;l ) Fig. 3: Bilinear form of the consolidation line with the definition of an instantaneous modulus i 3.1 Nurnerical solution using the collocation method A suitable method for solving Eq. (12) combines the collo- cation method alongwith a cubic spline approximation of the pore pressure distribution. After discretization and introduc- tion of instantaneous parameters i and i6 Eq. (12) becomes nf (t1a) = -x ^to.i(tj+r)rudty) + r uB;(t1*t)* (,a) +i1,(tp). Further suppose an equidistant spread ofthe collocation points with srep Ah = zi- zi_r. The cubic spline equation rhen reads l- r ? 2Mo-t * 2Mi + : M;*1= ft(0,-t -2F; + h*) . (19) Combining Eqs. (19) and (18) then readily provides the final tridiagonal system ofalgebraic equationsin the form lr ' I fi * s' a'"'- ft i a) ft lu,-'(rr*, ) * | . rl +f2 -6rArcrr\ttu) oor- lri(/r*,) * ,r0l [r r I * f , * s' rr " '"('1) fi lu +{t 1.t) = 8. = fi fr,-,t(, i) - 2PtQ j)+ Pi.,r(r; )], where pr(t ) = r u p;(t 1) * 11. (t pt) vo(t)= fi(rr(,,) - Fi-{t)) * * {Q*o(,)+ uo-,(t,)) To complete the numerical procedure the following boundary and initial conditions are supplemented in accor- dance with the experimental set up r{t1)=0, o,(t1)=0, (23) F;@)=-o-* o'{,010) i=2,...,n -r, where n is the number of collocation points. The selected cu- bic spline boundary conditions assume the form ( l5) rrrt(t)=M,(t,)=0. 3.2 Numerical solution using the finite aolume method It is evident that a direct solution of Eq. (12) Ieads, in gen- eral, to a nonlinear system ofequations. Therefore it appears preferable to start by discretizing the governing equations, which are then successively used in a single time step. In view of the one-dimensional problem, Eq. (9), the set of governing equations will be discretized using a three-layer model. The individual layers (finite control volumes) have different thick- ness, and diverse materials can be assigned to each layer. The pore pressure is assumed to be constant inside each control volume andwithin each time step, and to be equal to the pore pressure at the so called grid-point, Fig. 4. Its position is pre- scribed by a chosen parameter/. Trarufer eqrnti,on Application of the finite difference scheme converrs Eq. (4) for the water flux throughout a layer i with thickness Az. into the form (2 1) /99\ (14) Note that by using the above approximation we force Eq. (12) to be fulfilled only ar a cerrain points of colloca- tion (i); Y,(t) and M,(1,) then represenr the first and second derivatives -of the unknown pore pressure at the id point, respectively. AIso note that the pore pressure dependent material parameters K and e are found from an initial pre- diction of pore pressure T;'(tp) at time 6u r-q. i7). Accepting such assumptions allows the above equation to be written in the form bi (t 1.t) * "(t 1 a) rrr ;(ri*t ) = p(ri*, ) where o,(r7*,) =. at,#tP "f,(t i), F;(r;.r ) =. e' W* ;','l *rj f,o(, i)f , "fff(r7-r) = o^{i(tj)+ pyQjn)- fi(t,) To solve Eq. (15) we further assume the generalized trape- zoidal rule [5] and write ii (t p) = fi(of 1,,.r) - To' (t p)), T;'(t1-r)=Oi(t1,") + (r - )ttii(,i.t), where I = Q - t i) I u ; Lt = t i *1- I ; is the time step. parame- ter t should b6 thosen frorir r e(il2,l) to ensure numerical stability. The most common choice is r = f2. When substitut- ing Eqs. (16) and (17) into Eq. (15) we arrive ar l0 (94\ (16) (rry Acta Polytechnica Vol. 43 No. 112003 t,j |<-----+h------rh--------------- r L=i', liLz' (1-fi)Az, aZt*', Fig. 4: A three layer model ty =-+-(pi -pi-t)=-@.n\pi -FYt), eb)gAz;r" '- '' where | -r, . rgAzj-t , gLzi , =(r - JA)?+ ri _ (26) hzi K;_t K i is the overall resistance between grid points (i - l) and i. Balnnce eqm,tion Employing the generalized trapezoidal method provides the time variation and the rate of volumetric strain in the form, recall Eqs. (16) and (17), 1.t",;(t1*t)=;;(r,,'Qi.r) - 2,,,(,r.')), Q7) 2,,;(t1',) =r",i(ti)+(l - r) tt z,,i(t1) . (28) Next, applying again the finite difference scheme yields an incremental form of the balance equation (5) for layer I A. .nw , ff; (r".,(,r.')-2,,'(r;.r))* (zg) + 1ia(t 1 *r) - li' (t, *1) = o. State eqntion Introducing Eqs. (16) and (27) into Eq. (8) gives i-1 a i+1 I 3.3 Nurnerical solution using the finite element method In view of the general solution of fluid-solid interaction we recall the finite element method (FEM) often used to solve differential equations similar to Eq. (9). Apart from Eq. (9), which must be satisfied at any point inside the soil body, the boundary value problem requires formulation of the boundary conditions. To that end, the boundary is de- composed into two parts, | =fp +Tn, with the following boundary conditions: Essentinl boundary cond;ition p. = F., onl, , N atural boundary condition t:+€ * a!- * 4 =s, on rn , (s4) lu 'lw OZ where q represents the prescribed fluxes. Fig. 5 shows the boundary conditions pertinent to isotropic consolidation. In this particular case, no fluid is transferred across the interface f, and therefore g = 0 and also # =, p*= 0 onfp p*(t,z) q=0onfq Fig. 5: Boundary conditions (30) To derive a set of finite element equations we start by in- troducing the balance equation (9) and its natural boundary condition (34) into the principle of virtual work, written as if "+"7f ?(*+.l- i.lap,a,*Jl ^l,L oz( oz ) l (35) Integrating the first term in Eq. (35) by parts yields HH t 4!-! g- K';,ff ( a, -[ ap,, i"'a, * q6p"'(H )= 0. (36) J Az y"')" "' Az J00 The next step employs the usual approximation of the water pore pressure in the form bt=Nndn,r (37) N7, stores the element shape functions and d, is the of nodal pore water pressures. After substituting the lt (33) T 1., I I ,,,0(t1*r)- G,,i(r;*r) = #fr* ffiQrO1*t) -Ti'(t'-'))' After combining Eqs. (30), (29) and (25) we arrrve at a set of tridiagonal algebraic equations in the form t (hz i*)pi*1(t 1*r ) - | (h.i*r ) + (hz ;) + L - Lripo' . L(ti*tt'z) I rA r + os(thN z) ;6lrt1rr*r) + +1nz)pi-1(tp)= =^###ffi,*p,.,1 At the onset of consolidation the initial condition, similar to Eq. (23), reads -tT.-. P;(o)=-o, + ofi.!(o) i=\,...,n, (32) where a is the number of finite control volumes. (31) where vector Acta Polytechnica Vol. 43 No. l/2003 above approximation inro Eq. (36) we identify the following matrices, [6], pmnzability matrix H . = Iffi *iiu(N7)-rya,, *L,p =Y, 0 compressibtlily matnx U alr c =-lNfNodz. J 0 Note that the last term in Eq. (36) drops out due to the prescribed boundary conditions. Since Eq. (36) must be satis- fied for any admissible virtual pore pressure W' *. arrive at the following system of equations capQjn) *F(t1r1)ae(17*r) =0 , (40) gh,ere the material paramerers entering Eq. (38) to get marrix f (tr+r) are evaluated under the same assumptions as in Section 3.1. We usually require the material parameters, including the coeflicient of permeability, to be constant within a given ele- ment. Eqs. (7) and (10) suggest that the volumetric strain is constant as well. When assuming a linear approximation of the pore pressune, the vector ofnodal pore pressure values of the it element becomes ..7(o,o)'={apr,dpz},. (4r) The constant effective mean stress offr in the i'h element is then given by ot#=o_ @+d 4 Numerical results Before proceeding with the results derived fiom the pre- sented formulation we draw the readers attention to lig. 6, which manifests several drawbacks associated with the solu- tion of the consolidation problem when applied to rhe stan- dard Cam clay model. The solid line shows results derived ex- perimentally. The theoretical predictions were derived as- suming the material parameters listed in Thble l. The dash-dotted line follows from numerical calculations assum- ing the linear consolidation line (virgin soil) and the deforma- tion dependent coefficient of permeability, while the dashed line was found from the bilinear consolidation line, but keep- ing the coeflicient of permeability constant. Clearly, neither the deformation dependent variation of the coeffrcient of (3e) (38) (42) 6(L 5 I o th E o I For the sake of consistencv. the volumetric strain e-. and the water pore pressure p' should be of the same trde, (recall Eq. (8)). Here, the desirable consistency of the solu- tion is achieved by applying so-called selective integration to integrate the compressibility matrix (one-point Gauss' quadrature). Once the discretization in space has been performed Eq. (40) represents a set of ordinary differential equations in time, which can be integrated numerically, t2l. To proceed, con- sider again the generalized trapezoidal rule such that ioQp) = j;(ooQ,.,) -ap(tp)), (43) Vo(t1',) =ap(,i) +(r-r)ttio(t1) . (44) Substitution of Eq. (a3) into Eq. (40) finally leads to a sys- tem oflinear algebraic equations given by "rt" * "tF(r7.,)) dpQin)=]c 70Q1.). (4b) The initial condition prescribing the water pore pressure at the beginning of loading, I = 0 takes the form itap=d;p(O)=-cn+";fi(O), i=1,...,n, (46) where n is the number of elements. t2 Time [min] Fig. 6: Time variation of pore pressure Table i: Material parameters rc [kPa] l" [kPa] eg Ko lrrls] afitr g