COMPARISON AND APPLICATION OF DIFFERENT ELASTO-PLASTIC CONSTITUTIVE MODELS IN FEM ANALYSIS OF AN EXCAVATED SOIL SLOPE

A modified anisotropic elasto-plastic model is presented in this paper based on the frame work of Modified Cam-clay model and Original Sekiguchi-Ohta model by applying a new stress ratio parameter, through which the yield locus and hardening axis in the p-q plane were adjusted to give a more reasonable description of the practical excavation stress path from initial K0 state. Typical excavated soil slope was selected as the illustration example to carry out the finite element numerical analysis by applying four constitutive models (Original Cam-Clay, Modified Cam-Clay, Original Sekiguchi-Ohta and Modified Sekiguchi-Ohta). The calculated displacements of the slope with above different constitutive models were then compared with the measured infield data. It was found the field measured data agree better with the results calculated from modified model, which indicates that the modified anisotropic model is more suitable in the description of the stress path of excavated soils slope after K0 consolidation.


Introduction
The initial stress state of the natural deposited ground is anisotropic in the sense that the vertical stress is typically larger than the lateral stresses. Some construction activities, such as excavation or backfilling, induce anisotropy in the stress system. The behavior response of soils is affected significantly by stress-induced anisotropy (Yildiz et al. 2009;Atkinson et al. 1990). Neglecting the effects of anisotropy may lead to inaccurate predictions of soft clay response (Karstunen, Koskinen 2008;Oliveira et al. 2010;Taiebat, Dafalias 2010). Cam-clay model (Roscoe, Schofield 1963) is one of the most important models which were developed based on the critical state theory. The yield surface is bullet-shaped in the original Camclay model (CC) and it is an ellipse in the modified Camclay model (Roscoe, Burland 1968) (MCC). But in both models, the yield surface is aligned along the p axis and hardens isotropically. To simulate the soils' response from isotropic to anisotropic, various revised models have been proposed based on Cam-Clay model (Graham et al. 1983;Dafalias 1986Dafalias , 1987Korhonen, Lojander 1987;Thevanayagam, Chameau 1992;Newson, Davies 1996;Wheeler et al. 2003;Gnanendran et al. 2006;Cui, Delage 2010;Yin et al. 2010). Most of these models, however, do not take into account the combined effect of anisotropy and excavation stress path. Furthermore, the application of these models to practical geotechnical design usually needs more calculation parameters compared with Cam-Clay model, which require nonstandard labora-tory tests. Sekiguchi-Ohta model (Sekiguchi, Ohta 1977) (S-O) was developed based on critical state theory to simulate the stress and strain relationship of normally K 0 -consolidated clays, especially for excavation stress path. A non-negative normalized shearing stress parameter was applied to account for stress-induced anisotropy, by which the initial K 0 state is taken into consideration.
To overcome the shortage of MCC model (isotropically hardening along p axis) and S-O model (bullet shaped yield surface). A modified anisotropic elasto-plastic model is presented in this paper based on the frame work of Modified Cam-clay model and Sekiguchi-Ohta model by applying a new stress parameter, through which the yield locus and hardening axis in the p-q plane were adjusted to give a more reasonable description of the practical excavated stress path from initial K 0 state. Four elastoplastic models: (1) Original Cam-Clay (CC); (2) Modified Cam-Clay (MCC); (3) Original Sekiguchi-Ohta (S-O) and (4) Modified Sekiguchi-Ohta (MS-O), were incorporated into a finite element program. Typical soil slope crossing section of the East Route of the South-North Water Diversion Project, China, in Sandong Province was selected as the illustration example to carry out the finite element numerical analysis. The calculated displacements of the slope based on above four models were compared with the measured infield practical displacement data. The comparison showed that, the field measured data agree better with the displacement results calculated from Modified Sekiguchi-Ohta model (MS-O) model. It indicates that the modified anisotropic

Typical stress path of excavated slope
Naturally consolidated soil is essentially a crossanisotropic material (Arthur, Menzies 1972). The main reason for the anisotropy is that most natural soils have been subjected to one-dimensional consolidation under initial K 0 stress state. There are general two categories of loading modes in different construction process: backfilling and excavation, just as illustrated in Figs 1a and 1b.
According to the different locations of point A, B, C and D there will be different loading or unloading stress paths in both backfill and excavation process, as illustrated in Figs 1c and 1d. And in each path, K 0 stress state is the typical state for natural soil foundation. In the initial stress state before backfilling or excavation, the major principal stress is in vertical direction and the minor principal stress is in horizontal direction. Plane strain assumption is usually taken in analysis; the intermediate principal stress is in the direction of plane strain. During backfilling, loads are applied in the direction of major principal stress. However, the loading mode is much different in excavation process: loads are reduced in the direction of minor principal stress. The different loading modes of have remarkable influence on the stress and deformation analysis of slope. Especially for the excavation, this kind of unloading stress path is much different from the traditional loading stress path simulating in most triaxial testing. It is necessary to adopt a reasonable model to describe the practical excavated stress path (Zhang, Yin 2007).

Original Sekiguchi-Ohta anisotropic elastoplastical model
Many constitutive models have been developed to reflect the initial K 0 and excavation effects. Sekiguchi and Ohta (1977) proposed an anisotropic hardening elastoplastic model for clays, which is an extension of the original Cam-clay model to normally K 0 -consolidated clays, applying the bullet-shaped yield locus like the Original Cam-Clay model as Fig. 2 shows. This model has been widely applied in the finite element computation to predict the stress and deformation of earth structure in engineering practice in Japan. The Original Cam-Clay model is an elastoplastic soil model based on isotropic conditions. The yield locus for this model is defined using a bullet shape. The model assumes that the elastic shear strain is zero and the soil dissipates the applied energy by undergoing plastic shear strains. The main limitations of the Original Cam-Clay model is that the application of normality at the intersection of the yield curve with the The compression stress path after K 0 is along 1-4 and the extension stress path is along 1-2-3. Actually, a cross-anisotropic soil undergoing pure vertical compression (vertical major principal stress) would mobilize higher shear strength compared to that undergoing pure shear (major principal stress at 45°) or pure vertical expansion (horizontal major principal stress). This means that a cross-anisotropic soil will yield at a much lower value of deviatory stress in extension than that in compression (Oda 1981;Burland 1990). Many researchers' work have proved that the yield locus is no longer symmetrical about the mean normal effective stress (p) axis in p-q stress space (Lade, Nelson 1984;Dafalias 1986) and the stress path after K 0 consolidation, or as illustrated in Fig. 2, is not along 1-2-3, but along the dash line. An asymmetric yield curve implies that the stiffness and strength of a soil in the vertical direction is significantly different than that in the horizontal direction. It is important to give a clear description of the initial K 0 state influence in developing model. Especially, when undergoing excavation, the loading direction is different from usually applied major principal stress direction loading mode, which influences the magnitude of the mobilized shear strength.
Many researchers' work proved that the initial anisotropic consolidation will have a remarkable influence on the yielding locus shape in p-q space and lots of anisotropic models were developed (Oda, Nakayama 1989;Liang, Ma 1992;Callisto, Calabresi 1998;Li, Dafalias 2002;Ling et al. 2002). The Sekiguchi-Ohta model (1977) is an extension of original Cam-Clay model with consideration of the anisotropy by using a new stress ratio. The difference between the Sekiguchi-Ohta model and the original Cam-Clay model is only the stress ratio used in the in the yield function. In the Sekiguchi-Ohta model, the relative stress ratio * η , by which the initial anisotropic stress state is taken into account, is used instead of the stress ratio q p = η as in the Original Cam-Clay model, in order to express the rotation of the plastic potential. The yield function (f ) and the plastic potential function (g) for the Sekiguchi-Ohta model are expressed as follows: where: λ and κ are the slopes of normal consolidation and swelling lines, respectively; M is the slope of the critical state line; e 0 is the initial void ratio of soil for p = p 0 , p v ε is the plastic volumetric strain that is used as a hardening parameter, and the relative stress ratio * η is defined as: where: It is worth noticing that the material parameters ( λ , κ , M and v) in the Sekiguchi-Ohta model are the same as those in the Cam-Clay model. When shearing starts from the isotropic stress state: then there is: in which, 0 ij σ is the value of ij σ at the end of the anisotropic consolidation and ij δ is Kronecker's delta and: the Sekiguchi-Ohta model changes to the original Cam-Clay model, the yield function is:

Revised Sekiguchi-Ohta anisotropic elastoplastical model
Though S-O model is popularly applied in many fields, there are still some shortcomings existing: (1) the bullet shaped yield locus in p-q plane, which has the similar shortcoming as CC model. When simulating the isotropic consolidation initial stress state, the same problem occurs like the Original Cam-Clay model; (2) when unloading, hardening along K 0 line may induce smaller elastical zone and greater plastical strains compared with practice situation. The modified Cam-Clay model developed by Roscoe and Burland (1968) belongs to the class of critical state models formulated within the framework of classical plasticity theory. Because of its capability to account for the most relevant characteristics of soil, the MCC is one of the most widely used material models for geotechnical applications. In this model, it is assumed that the dissipation of energy is due to both the elastic and plastic shear strains. This assumption required the yield curve to be an ellipse. It is expressed by a yielding criterion represented by an ellipse oriented in line with the p axis and defined as a function of the pre-consolidated pressure, p 0 . The ellipsoidal yield surface in the p-q stress space give a good description of isotropic hardening and softening mechanism along p axis. To get a rotated harden axis and ellipsoidal yield surface, a new stress ratio '* η was introduced as: which means the yield locus hardens along the 0 / n η axis in p-q stress space, n could be 1, 2, 3 or other figures to represent the different angles between the harden axis and p axis. And for the model in this paper n was set as 2 to represent the harden axis from p axis to K 0 /2 line in p-q space, which gives a reasonable simulation of the initial stress state. And the yield function was written as: in which: For initial K 0 state, there are: Then: If the shear stress begin from isotropic state, K 0 =1, '* / q p = η , above yield function (12) changes to the yield function of Modified Cam-Clay model: The yield surface of MS-O and MCC are given in Figs 3a and 3b. Table 1 gives the comparison of four different models' yield function and harden axis.

Geometry and infield measurement
The South-North Water Diversion Project is a multi decade project to better utilize water resources available to China. The main idea of this project is to divert water from the Yangtze River in Southern China to Northern China. It was sub-divided into eastern, central, and western transport systems. The Eastern Transfer Route will divert water from the lower reaches of the Yangtze River in Jiangsu province to the Hai and Huai River Basins in Northern China. Yellow River water diversion project is a sub-project of Eastern Routine from the Dongping Lake, Shandong Province to the east line of the Yellow River water north of the section of trunk water supply works, which is the key to controlling the Eastern Transfer Route. A typical crossing section near K2+300 of the channel was selected to carry out FEM analysis in this paper. Geometry and soil layer distributions of this section are shown as Fig. 4.  (1 ) Ellipse shape, Hardening along 0 / 2 η axis

Fig. 4. Geometry and Soil layer distributions of slope K2+300 (m)
The mesh for the FEM calculation is shown in Fig. 5, in which the excavation process was divided into four successive steps. There are totally eight soil layers of the whole soil slope. Totally five load steps were simulated in the FEM calculation. The initial self-weight stress state was calculated at the first calculation step before excavation of soil slope. Four following steps were applied to simulate the actual channel slope excavation process, as simulated in the FEM mesh shown in Fig. 5, in which four layer meshes representing different designed excavation elevation in each excavation step. The slope was excavated from elevation of 40.50 m to 34.33 m in the four excavation steps. Field displacement measurement was performed during the excavation process. Six points of slope are selected to carry out the comparison of measured and calculated displacements. The node numbers of these six points in finite element mesh (which are 172, 176, 197, 216, 254 and 255) and according measurement location (ES2-2, ES2-1, ES2-3, ES2-4, ES2-5 and ES2-6) are noted in Fig. 6. Fig. 6. Instruments of in field Test (Section K2+200)

Soil parameters
The material parameters of each soil layer used in FEM analysis for the constitutive models are listed in Table 2. Soil parameters used in the analyses were determined based on a combination of the results of the subsurface investigation and laboratory testing.  The analysis program BCF was coded in FORTRAN language, which was developed in Hohai University by Prof. Yin for 2D finite element analysis on soil consolidation and stress-strain analysis. Different constitutive models could be employed conveniently in the FEM computer simulations to describe the stress and strain relationship of the soils. To simulate the practical channel excavation construction process, some special calculating technique was applied in the finite element analysis. The finite element meshes were constructed one time at the first calculation step, which is the initial state of the unexcavation, only initial gravity stress was calculated and there was no displacement. With the development of the excavation process, some meshes representing excavation parts were excluded and reverse node force were added to simulate the excavation caused force. At each excavation step, some parameters of those elements representing the excavation soil will be changed, such as that a small elastic modulus value (in this program, 10 kPa was taken) will be set to these elements. And, the reverse equivalent load to simulate the excavation unloading will be added at the according elements' nodes.

FEM results and comparison with field data
In Figs 7a to 7f the observed and computed vertical displacements at points 172, 176, 197, 216, 254 and 255, are plotted against loading step, from step 1 to step 5. In all the figures, "v_TEST" represents the observed results and "v_CC" (Original Cam-Clay model), "v_MCC" (Revised Cam-Clay model), "v_S-O" (Original Sekiguchi-Ohta model) and "v_MS-O" (Modified Sekiguchi-Ohta model) repent the results computed by the finite element method using the different constitutive models respectively. The positive displacement means expansion deformation (upward) and negative displacement means settlement deformation (downward). It can be seen from above figures that the vertical displacements computed using the Modified Sekiguchi-Ohta model agree with the observed results better that other models when there are settlement deformation, such as point 172, 176, 197, 216 and 254 illustrated in Fig. 7a to 7e, which are located in the excavated area. When there is negative vertical displacement at point 255 shown in Fig. 7f, original Cam-Clay and Modified Cam-Clay model give a better result compared with Sekiguchi-Ohta and Modified Sekiguchi-Ohta models because the stress path in this area is different from that in excavation area.
Figs 8a to 8d give the calculated vertical displacements contour in different excavated steps based on Modified Sekiguchi-Ohta model. From Fig. 8, it can be noted that with the development of excavation, the vertical displacements increase. The positive displacements distribute near the excavation area and the negative displacements distribute in the location where there is no excavation. There is some potential tendency that the whole slope behaviors like slide from the bank to the center. Thus necessary slope safety analysis has been carried out based on the work given in this paper. The corresponding work will be presented in other paper.

Conclusions
By applying the ellipsoidal yield shape and rotating the harden axis, a modified anisotropic elastoplastic constitutive model was presented in this paper based on the frame work of Original Sekiguchi-Ohta model and Modified Cam-clay model. By which the initial K 0 stress state and excavation stress path could be taken into account properly. Typical cross-section of Yellow River water diversion was selected as the illustration example to carry out the finite element numerical analysis. Four elastoplastical constitutive models were applied in the finite element numerical analysis and the calculated results were compared with the infield displacement measurement. There is better accordance between the vertical displacement calculated from modified Sekiguchi-Ohta model with the infield measured data, which proves that the suggested model is more suitable to describe soil stress and stain relationship under excavation stress path. Several conclusions could be summarized from this study: 1. A modified anisotropic elastoplastic model for clays was proposed, which is capable of describing deformation and strength characteristics of clays with the excavation stress path after initially anisotropic K 0 consolidation state.
2. The modified model parameters are the same in the Cam-clay model or the original Sekiguchi-Ohta model, which makes it convenient to carry out the comparison calculation.
3. By adjusting the stress ratio parametre * η , the hardening axis of proposed model in p-q plan can be easily rotated to any axis between the K 0 line and p axis in p-q stress space.
The methodology and numerical analysis presented provide a better understanding of stress state of practical layered soil excavated slope and can be used as a valuable reference for analysis of similar cases. It is also suggested that proper constitutive model should be applied in the FEM analysis, which should describe different stress paths of the excavation soil mass.