Numerical Impleme ntation of rsotropic Consolidation of Clayey Soils

This paper reports on implementation of several numerical techniques to solve a set of governing equations resulting from simple one dimensional isotropic consolidation of soils that behave according to the Cam clay model. Three different methods of solving the equations of consolidation, namely the collocation method, the finite volume method and the finite element method, are presented. Apart from evaluating their efficiency, which becomes particularly crucial when implementing these techniques in the framework of an optimization problem aimed at tuning the model parameters, a set of parameters of a Cam clay model driving the time dependent response of the soils (deformation dependent variation of the coefficient of permeability and preconsolidation pressure) is also discussed.


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 extensive 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 interaction.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 efficient, 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.Therefore, a critical evaluation of the efliciency of these methods, discussed in Section 3, is one of the goals of this paper.

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 -cF o*=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 Thansport of the liquid phase throughout the soil sample can be described by the following set of equations: 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.
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).
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 -.-, | , 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) 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 numerical 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 implemented 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 phenomenon, evident fiom Fig. 2 (see also [7]), results in a typical eo = 0.5, Ko = 2. 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 experimental measurements.
To avoid numerical difliculties around the point of discontinuity 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 initial prediction ofpore pressure taken at the end ofthe current time step.Details are given in Sections 3.1-3.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)  ( Combining Eqs. ( 19) and (18) then readily provides the final tridiagonal system ofalgebraic equationsin the form = fi fr,-,t(, i) -2PtQ j)+ Pi.,r(r; )], where pr(t ) = r u p;(t 1) * 11. (t pt) To complete the numerical procedure the following 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 thickness, 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 Note that by using the above approximation we force Eq. ( 12) to be fulfilled only ar a cerrain points of collocation (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).
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.

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 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 # =, To derive a set of finite element equations we start by introducing the balance equation (9) and its natural boundary condition (34) into the principle of virtual work, written as if "+"7f ?(*+.li.lap,a,*Jl ^l,L oz( oz ) l 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 J 00 The next step employs the usual approximation of the water pore pressure in the form bt=Nndn, r N7, stores the element shape functions and d, is the of nodal pore water pressures.After substituting the lt ,,,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 ;) +  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 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 solution is achieved by applying so-called selective integration to integrate the compressibility matrix (one-point Gauss' quadrature).
(44) Substitution of Eq. (a3) into Eq. ( 40 finally leads to a system oflinear algebraic equations given by "rt" * "tF(r7.,))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.In a nutshell, the first assumption essentially labels soils as being normally consolidated from their virgin state, which further contributes to a rather slow initial phase of the process of consolidation, Iig.6 (dash-dotted line).In their natural state, however, soils are usually in overconsolidated states, due to unloading from the original stress state (e.g., water table fluctuation).To account for the prior loading history we introduce the influence zone variation as a function of a certain level of the effective mean stress, usually referred to as structural strength.This level of the effective mean stress can be identified with preconsolidation pressure.In the present formulation, the structural strength is introduced through a bilinear shape of the consolidation line [8], [].The effect of this step becomes evident when examining Eq. ( 12).In particulal a lower value of l" at the initial stage increases the coefficient of permeability, leading to an acceleration of the consolidation process at this stage, a phenomenon observed in the experiments.
experiment --m=0,":f=-4okpa .-' m=0, aff=okea The second shortcoming comes fiom the originally as- sumed constant value of the coe{Iicient of permeability, re- sponsible for an abrupt increase in the rate ofconsolidation as the effective mean stress developed during consolidation becomes appreciable, Frg.6 (dashed line).Howeveq this stage of consolidation is associated with a significant deformation of the skeleton, leading to a decrease in the volume of the pores and subsequently to an decrease in the coeffrcient of permeability.Based on our extensive numerical investigations of isotropic consolidation [8], we proposed an exponential relation between the coeffrcient of permeability and the actual void ratio represented by Eq. (10).Note that this relationship was suggested by the authors based solely on numerical results.Nevertheless, it conforms with an experimental observation, [9] pp.138, fairly well.This set of numerical experiments also served to examine the individual numerical techniques.The results were rather inconclusive as to what method should be preferred, since all methods performed equally well in this very specific consolidation problem.Thus the fact that the finite volume method is usually more stable ruled the selection of this method for the subsequent numerical program, which provides a qualitative description of the above mentioned requirements.
First, the influence of the bilinear consolidation line trig.l, is investigated.The location of the slope discontinuiry on the rl, line corresponds to a given value of the structural strength of the soil.This value can be estimated using, e.9., the Casagrande method [3J.Tne results for various values of the structural strength offtare plotted in lig.7, assumrng tn equal to 6.It is evident that increasing the magnitude of ofitr speeds up the rate of consolidation in the hrst stage quite substantially.The second objective aims at slowing down the consolida- tion process in the second stage of consolidation associated with a decreasing volume of pores due to skeleton deforma- tion.This phenomenon is taken into account by incorporat- ing Eq. ( 10) into the numerical procedure.Results forvarious va$es of exponent m and for a given structural strength o#' = -30 kPa are presented in Fig. 8.The dash-dotted line shows a rather gradual decrease in the pore pressure profile with increasing time, which is in better agreement with experimental observations.Frg. 9 then summarizes the present theory.First, recall that the dotted line was obtained without taking into account the bilinear consolidation line.In such a case, a rather high value of rz locks up the numerical procedure at the beginning of consolidation.The dashed line then re ects incorporation of both requirements.Up to now a simple trial and error proce- dure was employed to derive optimal values of rz and of;" .Furthermore, reproducing laboratory data requires that we supply, apart fiom the structural strength parameters @li' ,* ), the following material parameters: the initial void ratio ag, the initial coeflicient of permeability Ko, the swell- ing index r and the compression index 1,.The last two material parameters and the starting value of precon- solidation pressure could be specified fiom the steady state response corresponding to the diagram shown in Fig. l.This, however, would require carrying out a number of isotropic consolidation tests for several predefrned levels of isotropic pressures attained at the end of the consolidation process, while measuring at the same time the change in volumetric strain Aer.Such an approach would become not only time consuming but also burdened by an additional experimental error associated with the difficulty of measuring Aer.Also, the determination of 7s and Kowould require additional labo- ratory tests.A simple solution to this problem, however, is offered by following the steps of mixed experimental and numerical methods.In such a case, combining the experi- inferring the desired parameters fiom a single laboratory test.This issue is discussed in [0].

Conclusions
The present paper was concerned with an effect of structural strength in the computational model of isotropic consolidation in which the skeleton deformation is governed by the modified Cam clay model.Several numerical rech- niques to solve the isotropic consolidation problem were explored.The finite volume method seems to perform best in the present problem of one dimensional isotropic consolidation.
Both the numerical results and the experimental data proved the need to incorporate the structural strength of soils and the time variation of the coefiicienr of permeability into the computational model.An independent experimental observation confirmed the exponential form of the relation- ship benveen the coefficient of permeability and the current void ratio proposed by the authors.
time derivative a( ) I at .

Fig. I :
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

Fig. 2 :
Fig. 2: Time variation of pore pressure with no smoothing of the bilinear consolidation line
Fig. 6: Time variation of pore pressureTable i: Material parameters

Fig 7 :
Fig 7: Time variation of pore pressure