BOUNDING SURFACE COUPLED FINITE ELEMENT CONSOLIDATION ANALYSIS OF NORMALLY AND OVERCONSOLIDATED CLAYS

The radial mapping version of the bounding surface plasticity model is implemented in a computer program to predict the response of cohesive soils. The eight-noded isoparametric element and Biot’s theory are used in this study for analyzing soil consolidation problems. The model has been used in the analysis for two classes of problems. The first involves the comparison of model predictions with the results of laboratory tests in compression and extension for normally and overconsolidated clays. The second class involves using the model to predict the results of one-and two-dimensional finite element problems of soil consolidation. The comparisons with experiments demonstrate that the model, through its simplicity, can describe realistically the soil response under different monotonic loading conditions at any overconsolidation ratio. The comparison between the bounding surface plasticity model with the classical modified Cam clay model shows considerably different rates and magnitudes of settlement


INTRODUCTION
Consolidation plays an important role in many soil mechanics problems.This is evident by the vast amount of literature devoted to the solution of this problem since the pioneering work of Terzaghi.Probably the most difficult problems that a soil engineer is asked to solve, is the accurate prediction of the settlement of a finite loaded foundation on a thick compressible stratum.
A major problem in applying any solution procedure to geotechnical engineering problems is to provide a realistic representation of the stress-strain characteristics for the porous medium.The choice of appropriate constitutive or stress-strain laws or models may have a significant influence on the numerical results obtained.Their importance has been enhanced significantly with the great increase in development and application of many modern computer-based techniques such as the finite element, finite difference, and boundary integral equation methods.It has been realized that the advances and sophistication in the solution techniques have far exceeded our knowledge of the behavior of materials defined by constitutive laws.As a consequence, very often, results from a numerical procedure that may have used less appropriate constitutive laws can be of limited or doubtful validity.
The objective of this study is to implement one of the several relatively new and very promising plasticity models for soils (the bounding surface plasticity model) in a general nonlinear analysis program to establish the accuracy and limitations of the model formulation.The model has been used in the analysis of two classes of problems.The first involves the comparison of model predictions to the results of laboratory tests.The second class involves using the model to predict the results of one-and two-dimensional finite element problems of soil consolidation.

THE GOVERNING EQUATIONS FOR SINGLE FLOW IN A DEFORMING POROUS MEDIUM
In this section, the governing equations for single flow in a deforming porous medium are developed.The governing equations are developed in line of Biot's self consistent theory (Biot, 1941;Biot, 1955;Biot, 1956).The solid phase is assumed to be comprised of a porous skeleton of particles surrounded by one fluid.The small strain theory is considered to be applicable, and so Darcy's law is assumed valid in terms of absolute fluid velocity.
It is assumed that a pure fluid pressure p causes only a uniform, volumetric strain by compressing the grains and that the major deformation of the porous skeleton is governed by the effective stress  .This is defined as follows, with the sign convention that tension is positive: (1) where is the total stress and m is equal to unity for the normal stress components and zero for the shear stress components.
The constitutive equation relating effective stress  to the strains of the skeleton is now independent of the pore pressure p, and for a general non-linear material can be written in a tangential form, thus allowing plasticity to be incorporated.If creep strain is present, the expression is written in a general form as (Lewis and Schrefler, 1987): where  d represents the total strain of the skeleton, represents the overall volumetric strains caused by uniform compression of the particles by the pressure of the pore fluid, with s K being the bulk modulus of the solid phase.
Finally, o  represents all other strains not directly associated with stress changes (swelling,  thermal, chemical, etc.).These types of strains are defined as "autogeneous" strains (Zienkiewicz et. al., 1977).
The matrix T D and the creep function c are dependent on the level of effective stress  and also, if strain effects are considered, on the total strain of the skeleton  (Lewis and Schrefler,  1987).
The equilibrium equation relating the total stress  to the body forces b and the boundary traction t ˆ specified at the boundary  of the domain  is formulated in terms of the unknown displacement vector u.Using the principle of virtual work (Zienkiewicz, 1977), the general equilibrium statement can be written as: for virtual displacement u  such that on the boundary part u  , where displacements are prescribed, these are not varied.Equation 3 is already a weak statement of the equilibrium relationship which also incorporates the boundary conditions.
The equilibrium statement (Eq.3) is also valid in incremental form: The effective stress relationship given by the equation 1 is now incorporated into this equation and the following expression is obtained: f d represents the change in external force due to boundary and body force loadings.
Further, on taking into account the constitutive relationship given by equation 2 and dividing by dt, the following equation is obtained: The geometrical complexity of a porous medium renders impossible a strict analytical treatment of the fluid velocity within the porous space.To overcome this obstacle, the fictitious seepage velocity (also known as bulk or Darcy's velocity) is defined as (Bear, 1972): where k is the absolute permeability matrix of the medium,  the dynamic viscosity of the fluid, p the fluid pressure,  the density, g the gravity, and h is the head above some arbitrary datum.The continuity of flow requires that the following expression is valid (Crichlow, 1977): which, on combining with Darcy's law given by the equation 8, results in: There are many factors which contribute to the rate of fluid accumulation and these are enumerated as follows (Lewis et. al., 1976): a.
Rate of change of total strain Rate of change of the grain volume due to pressure changes where n is the porosity.c.
Rate of change of saturation where s is the degree of saturation.d.
Rate of change of fluid density Finally, the change of grain size due to effective stress changes t    from equation 2a into 11e yields the following expression: The continuity equation for water (with no source term) therefore becomes: Equation 12 can be simplified by assuming water flowing at saturated conditions and dividing by  results in equation 11d for water being written as: where w K is the bulk modulus of water.Equation 12 then becomes:

FINITE ELEMENT APPLICATION IN SOIL CONSOLIDATION
The fully coupled solution of the one-phase flow equation in an elastoplastic porous medium will now be discussed in this section.The particular form of the continuity equation 14, together with the equilibrium equation 7, form the governing equations for soil mechanics problems within the line of Biot's self consistent theory.
The consolidation problem is a boundary value problem, and this type of problems requires that the governing equations are satisfied within all points of a continuum (domain  and that the boundary conditions are satisfied on the boundary of the domain.
The equilibrium equation 7 has the boundary condition already incorporated.Attention is therefore focused on the continuity equation.In this case the boundary conditions satisfy: a.
The continuity of flow across the boundary, where n is the unit normal vector and q is the outflow rate per unit area of the boundary surface.b.
Prescribed pore pressure.The condition that the continuity equation 14 applies throughout the continuum and that equation 15 applies on the boundary requires that: where a and b are a set of arbitrary functions since A and B are identically satisfied throughout their respective domains.
Conversely, if equation 16 is valid for any arbitrary values a and b, then the differential equations 14 and 15 must be satisfied at all points within and on the boundary of the continuum.
The finite element method will be applied to equations 7 and 16 in terms of displacements and pore pressures.In equation 14 the appearance of second derivatives for ) gh p (   necessitates a smooth distribution in space due to the integration of these variables.In order to overcome this limitation, a weak form of equation 14 is obtained by means of Green's theorem.
Upon substitution of equation 14 and equation 15, equation 16 becomes: Since the values of a and b are arbitrary, it can be made: a b   and thus eliminate some of the terms of the boundary integrals.Equation 17 therefore reduces to: The finite element approximation is now applied to equations 7 and 18.The displacements and pore pressures are expressed in terms of their values u and p at a finite number of points in space.
The expressions for u, p, and  take the form: in which N, N are the shape functions of displacement and pore pressure, respectively and B is the strain-displacement transformation matrix.Substituting equations 19 into equations 7 and 18, the finite element discretization gives the result: Equation 20 is valid for any value of the virtual displacement u  and can hence be written as: where in which K= the tangential spatial stiffness matrix; L = the coupling matrix representing the influence of pore pressure in force equilibrium; C = the creep matrix and df = load vector equivalent to the body force, surface traction and autogeneous strain, respectively.
The form of the function a in equation 21 is still quite arbitrary and must be specified before equation 21 can be solved.It is desirable to choose a form which will increase the accuracy of the approximation used.For this purpose, the method of weighted residuals procedure has been applied, using the Galerkin method (Zienkiewicz, 1977).
The function a is replaced by a finite number of functions within each element, which is in the Galerkin method are identical to the shape functions N .Equation 21 now becomes: where: in which H = is the spatial flow (or seepage) matrix; S = the compressibility matrix and f = the load vector equivalent to fluid flow of source elements, creep function and gravity load, respectively.
It can be easily verified that the complete set of equations is symmetric if the matrix T D is symmetric.
The integration of these equations usually requires the use of numerical techniques, and a standard method is that of Gaussian quadrature (Zienkiewicz, 1977), where the integrands are evaluated at specific points of the element and boundary surface and then weighted and assumed.The procedure is carried out in terms of a set of local coordinates  and  having values of 1  on the element boundaries.
Since the discretization in space has been carried out, equations 22 and 24 now represent a set of ordinary differential equations in time.For convenience, the equations are written in the following form: The values of u and p at different values in time may now be obtained by means of appropriate time-stepping algorithms.
The method used for the time discretization may be regarded as a one-dimensional finite element scheme as distinct from the spatial discretization, Kantorovich type approach, (Zienkiewicz, 1977).
The time domain is divided into a number of elements or steps and integration is carried out for each step to obtain the change of the parameters u and p .The step-by-step integrations may then be summed to determine the total change of the parameters.The integration takes the same form as used for the spatial integration, i.e., if F = 0 then   0 Fdt g where g is an arbitrary function of time.
When applying this method to equation 26, yields the following equations: where k t  is the length of the k-th time step.The first order-time derivatives of u and p may be approximated by assuming a linear variation of u and p within each step.
. The derivatives with respect to time of 1 N and 2 N are given as: By substituting equations 28 and 29, equations 27 take the form: Equation 30 may now be integrated, using the point collocation method, then divided by g , and finally rearranged in the form: Equations 31 are formed for all internal nodes of the domain and those boundary nodes where pore pressure value and/or displacement are not prescribed.The number of equations is thus equal to the number of unknown variables.
The complete set of equations may be used in the time-stepping procedure outlined above to determine the values of u and p at any point in time relative to their initial values.
The eight-noded isoparametric element is used in this study.The nodes in the elements used represent both displacements and pore fluid pressures, which vary quadratically over the element, as field variables on the same mesh.In other words, the nodes represent displacements and pore pressures at the same time as if there are two concurrent meshes.

NUMERICAL IMPLEMENTATION OF THE BOUNDING SURFACE PLASTICITY MODEL
Soil plasticity problems are nonlinear and history-dependent and thus require more elaborate solution schemes for boundary value problems than simple linear elasticity problems.All the techniques for nonlinear analysis can, with certain qualifications, be applied irrespective of the constitutive law, although some techniques are better suited to particular laws than others (Naylor and Pande, 1981).
The bounding surface plasticity model can be used to supply properties for most numerical solution schemes applicable to soil stress analysis problems.In general, the use of a history dependent constitutive model (such as the bounding surface plasticity) in a stress analysis program requires some form of an incremental solution procedure.In addition, in order to be able to employ reasonably sized solution steps (for implicit implementation), iteration within each step is also usually necessary.
The constitutive relation in inverse form is given by the following equation: The determination of the explicit form of the ijkl D tensor, equation 6 of reference (Dafalias,  1986) involves the use of equations 1, 7-10, 22, 28, 29, 35a, and 35b of reference.(Dafalias and Herrmann, 1986).Because the associated flow rule has been adopted, the ij R is equal to ij L and the I ' U is equal to I ' F , etc. Combining these expressions gives: In the above expressions, the elastic bulk modulus K is given by equation 29 of reference (Dafalias and Herrmann, 1986).The elastic shear modulus G is either defined independently or is computed from K and a specified value of Poisson's ratio.The plastic modulus p K is obtained from equation 30 of reference (Dafalias, 1986) in which the hardening function H ˆ and the bounding plastic modulus are defined by equations 30 of reference (Dafalias, 1986) and 28 of reference (Dafalis and Herrmann, 1986), respectively.) ( h is the heavy-side step function.The quantities related to the bounding surface such as b, I ' F , etc., are given in Appendix I of reference (Dafalias and Herrmann, 1986).The stress invariants are defined in equation 2 of reference (Dafalias and Herrmann, 1986).
The size of the bounding surface is controlled by the measure of the preconsolidation history o I .
The evolutionary law for o I is given by the second equation of 27b in reference (Dafalias and Herrmann, 1986), i.e.: The parameters and l I are described in reference (Dafalias and Herrmann, 1986).Using equation 26 of reference (Dafalias and Herrmann, 1986) and expressing the elastic volume change in terms of the bulk modulus K and the change in I gives: Dividing by o I and integrating over substep m gives: The first two integrals may be evaluated exactly, while the third is approximated by the trapezoidal rule (K is given by equation 29 of reference (Dafalias and Herrmann, 1986)).Integrating and solving for m o I yields: Integrating this expression yields:

EXAMPLES
In this section, in order to check the validity of the numerical implementation of the bounding surface plasticity model, a comparison is made with the experimental results of soft clay response under undrained monotonic deviatoric loading in compression and extension for normally and overconsolidated clays published by Banerjee and Stipho (1978;1979).
The solutions for one-and two-dimensional plane strain consolidation problems are shown in this section.In this case, a comparison has been made between the settlement predictions obtained by Siriwardane and Desai (1981) and those obtained by using the bounding surface plasticity model that is used in this study.The conventional modified Cam clay model was used by Siriwardane and Desai for the elastoplastic analysis.

Model Response and Comparison with Experiments
The predictions of the bounding surface plasticity model at the constitutive matrix level (Eq.33) and comparisons with experimental data for different overconsolidation ratios in compression and extension tests are now presented.
A series of comparisons between the experimental observations of the stress-strain in addition to pore water pressure and those predicted by the theoretical model are presented in Figs. 1, 2, 3, 4, and 5.These figures show by discrete symbols the experimental data for a triaxial undrained loading at overconsolidation ratios OCR = 2, 5, and 12 in compression tests and at OCR = 2 and 10 in extension tests and by continuous lines for the model simulation.All the theoretical results have been obtained by using the values of the input material parameters listed in Table 1.
The experimental data are taken as reported by Banerjee and Stipho (1978;1979) where samples of kaolin (LL = 52, PL = 26) were tested under undrained conditions.The method of preparation of these samples and details of the testing procedure are given by Banerjee and Stipho (1978).
Comparison with experiments demonstrates that the predicted response by using the bounding surface plasticity model is in good agreement with experimental data and that the model can describe realistically the soil response under different loading conditions at any overconsolidation ratio.

One-Dimensional Consolidation
In this section, the solution for a one-dimensional plane strain consolidation problem is shown.The results obtained from the bounding surface plasticity model are compared with those of Siriwardane and Desai (1981).Isoparametric eight-noded elements have been used instead of triangular elements that were used by Siriwardane and Desai.
The finite element mesh is shown in Fig. 6.The width of loading B is assumed to be equal to 0.102 m.An external surface load (p o = 47.9 kN/m 2 ) is applied at the top surface of the model.The input material parameters used in this problem are shown in Table 2. Some of these parameters are the classical material parameters within the critical-state soil mechanics context.So these material parameters are taken as reported by Siriwardane and Desai (1981).The others are the new material parameters related to the bounding surface plasticity model.The typical values of these parameters in the more limited ranges for practical applications are used in the analysis of this problem.Pore Water Pressure Fig. 7 shows the comparison between the calculated settlement at a typical node (107) using the bounding surface plasticity model with that from the problem that was solved by Siriwardane and Desai.The settlements are non-dimensionalized with respect to ultimate settlement predicted by Terzaghi's theory of consolidation.The non-dimensionalized time factor ( ) is used as an abscissa and E o is used to calculate c v .
The dissipation of pore water pressure at a typical node (31) and the vertical section (A-A) as shown in Fig. 6, are shown in Figs. 8 and 9, respectively.Here, the pore water pressures are nondimensionalized with respect to p o , which is the externally applied surface loading.
Comparison between the results from the bounding surface plasticity model with those of conventional modified Cam clay model by Siriwardane and Desai (1981) shows that the use of the first leads to higher settlements at initial time levels and smaller final settlements.Also, it can be seen that the use of the bounding surface plasticity model increases the magnitude of dissipation of pore water pressure especially in the final consolidation zone.Therefore, the predicted pore pressures are less than those from the conventional modified Cam clay model.

Two-Dimensional Consolidation
The solution for a two-dimensional plane strain consolidation problem is shown in this section.The same problem was solved by Siriwardane and Desai (1981).In this case, a comparison is made between the settlement predictions obtained by the bounding surface plasticity model and the results obtained by them.The critical state model was used by Siriwardane and Desai as the constitutive relationship.Triangular elements were also used to solve this problem while isoparameteric eight-noded elements are employed in this study.
The problem to be solved and the finite element mesh are shown in Fig. 10.The width of the loaded area B is assumed to be equal to 3.05 m.An external surface load (p o = 47.9 kN/m 2 ) is applied at the top surface as shown in the Figure .The problem is solved using the input material parameters shown in Table 3.The predicted pore pressures are non-dimensionalized with respect to p o , where p o is the externally applied surface loading.The initial modulus E o is used to compute the coefficient of consolidation c v value for use in the non-dimensional time factor T v .The predicted settlements are non-dimensionalized with respect to width of loading B.
Fig. 11 shows a comparison of timewise variation of surface settlement from the bounding surface plasticity model with those from a problem that was solved by Siriwardane and Desai (1981) using the conventional modified Cam clay model.It can be seen that the settlements from the two models do not differ significantly at initial time levels.However, at higher times the bounding surface plasticity model shows higher settlements but a smaller final settlement.Dissipations of pore water pressure at section (B-B) in Fig. 10 are compared in Fig. 12.The bounding surface plasticity model shows higher dissipation of pore water pressures at earlier times, while, at higher time levels, the two results tend to be similar.

Table 1 :
The Input Parameters Used to Predict the Response of the Bounding Surface Plasticity Model

Fig. 6 :
Fig. 6: Finite Element Mesh for the One-Dimensional Consolidation Problem

Fig. 9 :Fig. 11 :Fig. 12 :
Fig. 9: Pore Pressure versus Depth along Section A-A at Three Different Time Values Using the Bounding Surface Plasticity Model 059 Bounding Surface Plasticity Model Tv = 0.059 Modified Cam Clay Model Tv = 0.504 Bounding Surface Plasticity Model Tv = 0.504 Modified Cam Clay Model Tv = 1.06 Bounding Surface Plasticity Model Tv = 1.06 Modified Cam Clay Model

Table 2 :
Input Material Parameters for the One-Dimensional Consolidation Problem