- Email: [email protected]
Multi-stage finite element algorithm for excavation in elastoplastic soils
Compurers & Slrucrures Vol.46,No. 2, pp.289-298, 1993 Printed in GreatBritain.
004s7949/93 $6.00 + 0.00 0 1993 Pergamon Press Ltd
MULTI-STAGE FINITE ELEMENT ALGORITHM FOR EXCAVATION IN ELASTOPLASTIC SOILS E.
COMODROMOS,
TH. HATZIGOGOS and K. PITILAKIS
Department of Civil Engineering, Division of Geotechnical Engineering, Aristotle University of Thessaloniki, Thessaloniki, Greece (Received 12 November
1991)
Abstract-A
multi-stage finite element algorithm for simulating excavation in elastoplastic soils is presented. Based on the conception of the variable domain vector, the non-linear finite element equations are derived in order to take into account the appropriate domain of interest only. Numerical instabilities arising from element and node subtraction during the numerical procedure are overcome by using the static condensation algorithm in the equation solution. Numerical examples demonstrate that the proposed method satisfies the basic criteria of stress-free excavation surface and the principle of uniqueness for linear elastic materials. Furthermore, the notion of the variable domain vector in conjunction with the static condensation algorithm leads to the elimination of any contribution of the removed elements in the numerical procedure and consequently guarantees no violation of the principle of virtual work resulting from the simulation process. In the case of an elastoplastic analysis, a dependency of the provided final solution upon the expanse of the plastified zone has been found. In that case, the one stage analysis is unable to take into account modification arising from topological variations of the structure, such as tangential elastoplastic constitutive matrix. A final solution closer to the exact one is achieved by a multi-stage excavation due to intermediate stages which provides a more accurate behaviour of the previously plastified element.
INTRODUCTION
zero nodal force vector acting on removed nodes, proposed a hybrid method for calculating nodal forces which may preserve uniqueness, but at the expense of added analytical complexity. Borja et al. [13] proposed a variational formulation-based method which in conjunction with the conception of an infinitesimal stiffness satisfies the principle of uniqueness. Finally, Comodromos et al. [14-l 61, have proposed a particular algorithm for nullifying nodal forces arising from excavated elements, which in conjunction with a double pivoting algorithm for the equilibrium solution process satisfies both the principle of uniqueness and the stress-free surface. The simplicity of the conception of the infinitesimal stiffness matrix for the removed elements used by the aforementioned algorithms is offset by the two following disadvantages:
The conception of the early numerical methods for simulating excavations [l-5], was based on the physically coherent and widely adopted hypothesis of a stress-free surface resulting from the governing equations and the boundary conditions of the equilibrium of a continuum body.
An additional criterion of efficiency, accuracy and acceptance known as the principle of uniqueness was postulated by Ishihara [6]. According to the principle of uniqueness, which arises from the principle of superposition, there exists a unique solution independent of the sequence of excavation for a linear, time-independent and elastic material. Ghaboussi and Pecknold [7] demonstrated the validity of the principle by employing a virtual work formulation over the appropriate domain of interest. The theoretically ideal uniqueness is not achieved by many widely used excavation algorithms such as those proposed by Christian and Wong [8], Clough and Mana [9], and Mana [lo, 111. On the contrary, numerical examples of the aforementioned algorithms outlined a significant dependency on the number of excavation stages and discretization of the problem domain. The main conception of the above algorithms consists of limiting the contribution of the excavated elements to the solution procedure by considering an infinitesimal stiffness for the removed elements and trying to nullify the forces acting on removed nodes. Desai and Sargand [12], to overcome the violation of the principle of uniqueness resulting from a non-
1. Violation of the principle of virtual work arising from the contribution of the excavated elements to the equilibrium solution process. 2. The whole procedure leads to a global stiffness matrix presenting infinitesimal pivots corresponding to removed nodes. The use of infinitesimal pivots produces a considerable numerical error which may be overcome by a particular double pivoting algorithm, but at the expense of doubling the solution time. The objective of this paper is to propose an efficient, finite element based, numerical procedure for simulating multi-stage excavation for soils 289
obeying an elastoplastic law with particular focus on the satisfaction against the violation of the principle of virtual work resulting from the simulation of the excavated region. GENERAL
FORMULA~ON
Consider the equilibrium of a body given in Fig, 1, The external forces acting onto the body are surface traction f” and body forces P. The equilibrium solution must satisfy the boundary conditions given by eqns (2) and (3) while at every point within the body eqn (1) must be satisfied too
The presence of a given finite element, e, in the domain of interest, for a given stage, can be arranged by the use of the notion of the variable domain vector (VDV) @,, . The VDV is a vector of scalars which in conjunction with the assemblage process, eqns (5) and (6), determine in a unique way the appropriate domain of interest at a given stage n
~2, = i r=l
Vu(u) -P u=ug
= 0
in f.!
on r,
a(h) = P on f,,
(1) (2) (3)
where V denotes the vector gradient operator, u is the Cauchy stress tensor, up is the prescribed displacement vector, Q is the problem domain, and T,uT, = r are the problem boundaries. In the context of the finite element method, the original problem domain arises from the assemblage of the m finite elements in which the domain has been subdivided. The appropriate domain for a given stage can be arisen from an appropriate assemblage of finite elements. Let e be a finite element of the body given in Fig. 1
12”,= f c= I
[’
dR’ = T@,$, Jn.a
(6)
in which ‘es,, = {fI,&, . . . , S&, . . . ,8&) is transpose of the variational scalar vector with 8:, = 1 when a >n otherwise B:, =0 (Fig. 1) and, Sz, = {Q’, . . . ,ff’, . . . , fy”) is the finite element domain vector. The principle of virtual work is affected by the aforementioned variations and becomes l+z, = YXR>
(7)
where the internal virtual strain energy is given by eqn (g),
The substitution
of eqn (5) to eqn (8) leads to
(4) When applied to excavation problems, both boundaries and domain of interest are stage dependent and eqns (l)-(3) are generalized to allow their variations with excavation stages while eqn (4) is modified in order to take into account the appropriate domain corresponding to a given stage n. The variation of the domain of interest, resulting from the excavation process, involves the subtraction of finite elements according to the stage of the excavation.
Area to be excavated
in which TV’,u’ and @&are the transpose of the strain vector, the stress vector and the variable scalar of element e, respectively. Within the framework of the finite element method, the element strains and displacements for a multi-stage non-linear analysis can be expressed as functions of the global displacements vector ‘II; = N’UI:
(10)
‘6; = B’UZ,
(11)
where N’ denotes the shape function matrix of element e, B’ the stress-strain transformation matrix and Ut the global displacement vector corresponding to stage n and iteration solution k. Furthermore, using the notion of the initial stress method which is convenient for an elastic-plastic analysis 117, 181, the final stresses at kth iteration of stage n, can be expressed as functions of the previous stresses and the actual strains by means of the following equation Fig. I. Schematic variation of domain and boundaries produced by an excavation.
Multi-stage
FE algorithm for excavation in elastoplastic soils which can be rewritten incremental form
in which
ec:=g
the well-known
Kkuk=R” ” ” ex,-R!’ IfI,cr” .
is the consistent tangential matrix obtained by evaluating the variation of the stress vector ‘us with respect to the strain vector ‘.sk If we substitute ind’eqn (9) for the element displacements, strains and stresses, using eqns (IO)-( 12) we obtain
e;,,‘u;
under
(13)
n
Win,= i
291
(18)
The matrix Ki is the global tangential stiffness matrix arising from the assemblage of the tangential element stiffness matrix ki O;,aTBe‘Cy:Be dR’,
(19)
where I” is the residual force vector resulting from the subtraction of the internal force vector R”,, from the external force vector given by eqns (20) and (21), respectively
Been, _ , dR’
e-l s C%
It is obvious that eqn (14) leads to a zero internal virtual strain energy for every element presenting n > a, i.e. when the element is subtracted. The external virtual work for a continuum body with no variations of the domain of interest is given by the well-known equation
+f
0 E,,‘Be ‘a, _ , da’.
R”,,= 2 e=l
TNefSdT’ + F’ (20) e=l I P (21)
Im
NUMERICALIMPLEMENTATION
Applying the notion of the VDV, eqns (5) and (6), and using for the instance F’ and r’ instead of Fj, and r; eqn (15) becomes tl;,, “U; TNefbdR’ P=l s RI
W,, = i
rU;WTNef”dI-e+ ZT;F’.
(16)
It can be seen, that in the case of a subtracted element, the first term of eqn (16) is equal to zero, while the last two are not; that violates the principle of virtual work. However, as will be mentioned later, the component of the external force vector in which the contribution of those terms is stored is eliminated by the static condensation solution process. Finally substituting eqns (14) and (16) to eqn (7) and multiplying from the left by the transpose of the global displacement vector, ZTi, we get the following governing equation
TNef=dT’
(17)
By directly invoking the variable domain vector in the variational formulation the appropriate domain of interest can automatically be obtained without having to renumber nodes and elements of the original mesh. From eqn (19) arises that ‘k, = 0 if a 2 n; the contribution of the initial stress vector became zero in that case as well. However, an equilibrium violation results from the surface forces vector and the concentrated loads of the removed elements. This violation and the mathematical problems arising from the semi-definite form of the global stiffness matrix due to the stiffness matrices of the removed elements can be overcome by the use of the static condensation algorithm in the solution of the equilibrium equations. Static condensation [ 19,201, is employed to reduce the number of element degrees of freedom corresponding to the removed nodes and thus, in effect, to perform part of the solution of the total finite element equilibrium equation system. For each excavation stage an assemblage of the global stiffness matrix and forces vector takes place based on eqn (17) and then a reformation reduction is performed according to the following process: if K,, = 0, then the row and column i of the global stiffness matrix and the row i of the force vector are eliminated as well. The outlined procedure leads to a subdivision of the problem in two independent substructures. The first one concerns the remaining domain and the second one the removed region. The variational equations for the remaining domain are not affected
292
E.
~~ODROM~
by the removed domain and so the problem can be solved with respect to governing equations and modified boundary conditions without any violation. The proposed method has been implemented in a FORTRAN program called nonlinear finite element analysis in geomechanics (NFEAG). It has been verified by comparing predictions with closed-form solutions and other numerical schemes for problems in simulation of excavation sequences. In order to demonstrate the validity of the whole numerical procedure two typical numerical examples are presented. The first, concerns a numerical patch test on a linear elastic one-dimensional excavation problem solved by Clough and Mana [9] and Desai and Sargand 1121. The second one, concerns a twodimensional plane strain excavation problem which has been solved for linear elastic behaviour and an elastic-plastic behaviour with or without hardening. Example 1. Simulation of a linear elastic one-dimensional excavation In order to verify the validity of their algorithm, Clough and Mana [9, 111, performed a linear elastic one-dimensional excavation (Fig. 2) and compared their results with those of the exact solution given below. They removed the upper 20 m by using a oneand eight-element simulation, using eight-node elements. The elastic parameters have been taken as follows: E = lO,O~ t/m2; y = 1.00 t/m2; Y = 0.2 and & = 0.5; the Young’s modulus, the unit weight, the Poisson’s ratio and the at-rest earth pressure, respectively. The same problem was solved by the hybrid method [12], with one-element excavation mesh. The exact solution can be provided by the following linear elastic stress-strain equation
in which f,, and rIij denote the volumet~~ strain and the Kronecker’s delta, respectively.
Area to be excavated
t
A
A
a
Fig. 2. One-dimensional excavation problem.
et al.
Fig. 3. One- and four-stage simulation of the one-dimensional excavation problem.
Taking into account the boundary condition and simulating the excavation by a surface traction of 20 t/m, eqn (22) becomes a
E “=~+(l+v)(l-2v)~~~
vE (23)
from which results a vertical displacement of 3.60001 cm. The same problem was simulated, using the method proposed here, by one and four excavation stages (Fig. 3) by removing one row-element at a time. The solutions arising were identical providing a displacement of 3.60001 cm (exact solution) for the nodes lying on the excavation surface. The vertical displa~ment profiles of the above solutions are given in Fig. 4. It can be seen that the algorithm proposed by Clough and Mana produced a solution highly dependent upon the number of excavation elements. On the contrary, the hybrid method is not affected by the number of excavation elements and stages and provides a solution very close to the exact one. However, both algorithms violate the principles of the one-dimensional form of the problem providing different displacements for nodes lying on the excavation surface. The proposed method, by taking into account the approp~ate domain of interest, eliminates any numerical parasitic error resulting from the cont~bution of the excavated elements and consequently provides the exact solution without any numerical error. Example 2. Simulation of a two-dimensional plane strain excavation The aim of this example is to verify the effeo tiveness of the proposed algorithm in the case of an elastoplastic analysis and furthermore to draw some useful conclusions resulting from such a case. Our interest is focused on an open vertical excavation given in Fig. 5. The mesh is composed of 72 eight-node isoparametric elements and a condition of plane strain is assumed. Nodes along the bottom are fixed in both directions, while nodes along the right-hand and left-hand boundaries may translate freely along the boundaries but are fixed against displacements normal to these
Multi-stage FE algorithm for excavation in elastoplastic soils
293
Drucker-Prager Failure Surface
\< \
/ k
i/ \
/
‘\ \.~
/ ,.:
l. \ ‘.
Fig. 6. Schematic presentation of the cone cap model.
i
Table 1. Summary of parameters used in second example
,,.,‘oneElemcnl . .., . .._+._.A’
Mlnr
1111
Constitutive law Y
Fig. 4. Example 1, vertical displacement excavation surface.
profiles of the
boundaries. In the analyses performed the elements in the upper right-hand corner, bounded by the thick line, were removed in one step and three steps, three row-elements at a time. The aforementioned excavation scheme has been analysed for linear elastic behaviour, perfect elastic-plastic behaviour (Drucker-Prager failure surface), and elastic-plastic behaviour with hardening (cone cap model) (Fig. 6). The parameters used in various analyses, given in Table 1, correspond to the Boston Blue clay [21]. W, D and R are parameters of the elliptical hardening surface, a and k are the Drucker-Prager failure surface parameters and K, G are the bulk and shear moduli, respectively. For elastoplastic analyses the convergence is achieved once the following criteria are satisfied:
Fig. 5. Two-dimensional
@N/m’)
rc,
K (kPa) G (kPa) k (iPa) W D (kPa-‘) R
Linear elastic
DruckerPrager
19.8 0.9 4700 2200 -
19.8 0.9 4700 2200 0.25 10
-
-
CAP 19.8 0.9 4700 2200 0.25 10 0.3 6.1 x 1O-4 4.13
1. Displacement criterion
(24) 2. Residual force criterion
llr2II - II4II < L llrl II ’ f’
(25)
where cdand q/denote the error tolerances which were taken as 10m3 and 10e5, respectively and [I./I is the Euclidean norm.
plane strain open excavation mesh.
294
E.
COMODROMOSet
al.
Linear elastic analysis
Single-stage and three-stage excavation analyses were performed. In the first, the elements 54, 55, 56, 62, 63, 64, 70, 71 and 72 were removed simultaneously while in the second one, a three-step layerby-layer excavation was considered. Figure 7 shows the evolution of the horizontal displacement profile at the vertical boundary of the excavation, line AB on Fig. 5, for the three-stage excavation analysis (configuration R, -+ R, -+ R, + Q,) and the profile corresponding to the single-stage analysis as well (% +RJE). It can be seen that the resulting final solution for both single- and multi-stage analysis are identical, a fact that verifies the principle of uniqueness. Drucker-Prager
= al, + &-k
= 0,
(26)
where 1, denotes the first invariant of stress tensor, J, the second invariant of deviatoric stress tensor while a and k are material constants related to the frictional and cohesive strengths of the material. Strains are elastic for stress paths remaining under the ultimate yield surface, but are elastoplastic if a stress path violates the surface. The same analyses as in linear elastic case were performed; the evolution of the horizontal displacement at the vertical excavation boundary is given in Fig. 8. A small difference between the two final solutions arises from the profile comparison. The
Steo Level
First
Second
Step Level
Third Step Level c
50 mm -
Step Level
Third SteD Level
-
criterion
In this case, an elastic-ideal plastic soil behaviour is considered according to the Drucker-Prager criterion given by the following equation f(a)
First
50
mm
-8-
First
-s
Second
+
Third
Step
+
lnltlsl
PosItIon
-
Step Step
Fig. 8. Example 2, evolution of horizontal induced displacements along the vertical excavation boundary produced by the Drucker-Prager analysis.
maximum difference, observed at node 173, is 1.5% of the displacement corresponding to R, configuration. The plastified integration points at the end of the excavation are given in Fig. 9. Cone cap model
Cap models describe the yielding behaviour of soil with an ultimate yield surface that is fitted with a movable end cap. Both the ultimate yield and cap surfaces are symmetric about the hydrostatic axis. As ultimate failure surface of the cone cap model, the Drucker-Prager surface is used while the strain-hardening elliptical end cap yield surface is described by eqn (27). The movement of the cap is controlled by the hardening function given by eqn (28). Strains are elastic for stress paths remaining within the region bounded by the ultimate yield and cap surface, but are elastic-plastic if a stress path violates any of the two surfaces. The model presented in Fig. 6 consists of a cone shaped Drucker-Prager ultimate surface, fitted with an elliptical end cap, eqns (26) and (27), respectively
./-(I,, ,/&, K) = 111-
L(K)]* + R*J: + [x(K) - L(K)]2
+ +
Flmt
-&
Sacond
(27)
Single-stage K =
Step
& = w[l - eXp(-DX(K)],
(28)
Step
+
Third
Step
+
Inltlal
POIltlon
/
Fig. 7. Example 2, evolution of horizontal induced displacements along the vertical excavation boundary produced by the linear elastic analysis.
where R is the ratio of the major to the minor axis of the elliptic cap, X(K) and L(K) define the intersections of the cap with the I, axis and the failure surface, respectively and K is the hardening parameter, & denotes plastic volumetric strains and W, D are material constants [21,22].
295
Multi-stage FE algorithm for excavation in elastoplastic soils
ii
I I
I
I
I
Fig. 9. Example 2, state of stress at integration points produced by the Drucker-Prager
The horizontal displacement profiles at the vertical excavation boundary corresponding to single- and multi-stage excavation analyses are given in Fig. 10, while the evolution of the plastification zone corresponding to Q,, sl, and i& con~gurations are given in Figs 1l-13. The difference between horizontal displacement at node 173 provided by single- and multi-stage excavation is 2.5% of the total displacement. Remarks
A comparison between the aforementioned displacement profiles for linear elastic analysis, the Drucker-Prager elastic-ideal plastic analyses and the cone cap model is presented in Pig. 14. Figures 15 and 16 show the displacement history of nodes 168, 173 obtained from the single-stage excavation (Q,, --) C&) multi-stage and the three-step excavation (Q, -+ a, -+ a2 -+ 12,). It can be seen that the cone cap model gives the greater displacement and the greater difference between single- and multi-stage analyses. This is due to the expense of the plastified zone arising from the delimitation of the linear elastic zone by the elliptical end cap. The accuracy and the capability of the proposed method for solving non-linear excavation problems with particular application to elastoplastic soils with or without hardening outline its effectiveness and usefulness.
instabilities arising from element and node subtraction during numerical procedure are overcome by the use of the static condensation algorithm in the equation solution. Numerical examples demonstrated that the proposed method satisfies the basic criteria of stress-free excavation surface and the principle of uniqueness. Furthermore, the notion of variable domain vector in conjunction with the static condensation algorithm eliminates any contribution of the removed eiements in the numerical procedure and consequently guarantees no violation of the principle of virtual work resulting from the simulation process.
A multi-stage finite element algorithm for simulating excavation in elastoplastic soils has been presented. Based on the conception of the variable domain vector, the non-linear finite element equations are derived in order to take into account only the appropriate domain of interest. Numerical
Step Level
First
bird Step Level -
+ CONCLUSIONS
analysis.
50 mm -
Singlestage Step
-8-
Flrllt
+-
Second
-g
Third
Step
-4%
lnltlsl
POIItlon
Step
Fig. 10. Example 2, evolution of horizontal induced displacements along the vertical excavation boundary produced by the cone cap analysis.
gelastic mc,
~Cotncar
~Tenslm
[‘ZiaFailure
Fig. 11. Example 2, state of stress at integration points produced by the cone cap analysis at stage one.
Fig. 12. Example 2, state of stress at integration points produced by the cone cap analysis at stage two.
Fig. 13. Example 2, state of stress at integration points produced by the cone cap analysis at stage three. 296
Multi-stage FE algorithm for excavation in elastoplastic soils
297
modification arising from topological variations of the structure, such as tangential elastoplastic consti-
tutive matrix. A final solution closer to the exact one is achieved by a multi-stage excavation due to intermediate stages which provides a more accurate behaviour of the previously plastified element.
50
-yt
LItmar-Elaatlc
+-
Druckw-Prager
-a-
conm Cap
-+
lnltlsl
REFERENCES
mm
Pomltlon
Fig. 14. Example 2, comparison between the linear elastic, the Drucker-Prager and the cone cap final solutions.
0
1
L
3
4
Horizontal
5
6
Dfsplacement
7
6
(cm)
Fig. 15. Comparison between the displacements of the node 168 produced by the linear elastic, the Drucker-Prager and the cone cap analysis.
-1
4
0
1
tioYizon:ai D:splafxmk7t
ikn)
0
9
IO
Fig. 16. Comparison between displacements of the node 173 produced by the linear elastic, the Drucker-Prager and the cone cap analysis.
In the case of an elastoplastic analysis, a dependency of the provided final solution upon the expanse of the plastified zone has been found. In that case, the one-stage analysis is unable to take into account CAS 46124
1. J. M. Duncan and C. Y. Chang, Nonlinear analysis of stress and strain in soils. J. Soil Mech. Found. Div., ASCE 96, 1629-1653 (1970). 2. J. M. Duncan and P. Dunlop, Slopes in stiff fissured clay and shales. J. Soil Mech. Found. Div., ASCE 95, 81-104 (1969). 3. G. W. Clough and J. M. Duncan, Finite element analyses of Port Allen and Old River locks. University of California, Berkeley, Report TE-69-3 (1969). 4. G. W. Clough and J. M. Duncan, Finite element analyses of retaining wall behavior. J. Soil Mech. Found. Div., ASCE 97, No. SM12 (1971). 5. P. Dunlop and J. M. Duncan, Development of failure around excavated slopes. J. Soil Mech. Found. Div., ASCE 96, 471-493 (1970). 6. K. lshihara, Relations between process of cutting and uniqueness of solutions. Soils Foundations 10, SO-65 (1970). 7. J. Ghaboussi and D. A. Pecknold, Incremental finite element analysis of geometrically altered structures. Int. J. Numer. Meth. Engng 20, 2051-2064 (1984). 8. J. T. Christian and I. H. Wong, Errors in simulation of excavation in elastic media by finite elements. Soil Foundations 13, l-10 (1973). 9. G. W. Clough and A. I. Mana, Lessons learned in finite element analysis of temporary excavations. In Proc. Zhd Int. Conf on Numer. Meth. in Geomech., ASCE, Blacksburg, VA (1976). 10. A. I. Mana and C. W. Clough, Prediction of movements for braced cuts in clay. J. Geotech. Engng Div., ASCE 107, 759-777 (1981). 11. A. I. Mana, Finite element analysis of deep excavation behavior in soft clay. Ph.D. thesis, Stanford University, Stanford, CA (1978). 12. C. S. Desai and S. Sargand, Hydrid FE procedure for soil-structure interaction. J. Geotech. Engng, ASCE 110, 473-486 (1984). 13. R. I. Borja, S. R. Lee and R. B. Seed, Numerical simulation of excavation in elastoplastic soils. Inf. J. Numer. Anal. Meth. in Geomech. 13, 231-249 (1989). 14. E. Comodromos, T. Hatzigogos and K. Pitilakis, Finite element excavation in elastoplastic soils. In Proc. &d Eur. Spec. Conf: Numer. Meth. Geotech. Engng, pp. 573-584. Santander (1990). 15. E. Comodromos, A finite element algorithm for simulating excavations in elastoplastic soils. 5rh Young Geotech. Eng. Co& Grenoble (1991). 16. E. Comodromos, Contribution to the analysis of excavation by the finite element method. Ph.D. thesis, Aristotle University of Thessaloniki (1991). 17. 0. C. Zienkiewicz, S. Valliapan and I. P. King, Elastoplastic solutions of engineering problems-initial stress finite element method. Int. J. Numer. Meth. Engng 1, No. 1 (1969). 18. C. S. Desai and F. J. Abel, Introduction to the Finite Element Method. A Numerical Method for Engineering Analysis. Van Nostrand Reinhold, New York (1972). 19. K. J. Bathe and E. L. Wilson, Numerical Methoak in Finite Element Analysis. Prentice-Hall, Englewood Cliffs, NJ (1976). 20. E. L. Wilson, The static condensation algorithm. Int. J. Numer. Meth. Engng 8, 199-203 (1974).
298
E. CO~OD~O~~ et al.
21. D. N. Humphrey, Design of reinforced embankments.
Ph.D. Thesis. School of Civil Engineerinn, Perdue University, West Lafayette, IN (1986). -. 22. W. F. Chen and G. Y. Baladi, Soil Plasticity-Theory and Implementation. Elsevier, New York (1986). 23. C. Y. Chang and J. M. Duncan, Analysis of soil movement around a deep excavation. J. Soil Mech. Found. Diu., ASCE 96, 1655-1681 (1970). 24. V. S. Chandrasekaran and J. W. King, Simulation of excavation using finite elements. Tech. Note, J. Soil Mech. Found. Div., AXE 100, 1086-1089 (1974).
25. C. S. Desai, Soil-structure interaction and simulation problems. In Finite Element in Geomechanics (Edited by G. Gudehus), pp. 209-250. John Wiley (1977).
26. C. S. Desai and H. J. Siriwardane, A concept of correction functions to account for non-associative characteristics of geologic media. Int. J. Numer. Anal. Meth. Geomech. 4, 377-387 (1980).
27. G. Dhatt and G. Touzot, In Une Prbentation de la MPthode des Uments Finis, Znd Edn (Edited by S. A. Maloine). Collection de I’Univercit~ de Compiegne, Paris (1984). 28. D. F. Stolle, An interpretation of initial stress and strain methods, and numerical stability. Int. J. Numer. Meth. Engng 15,399-416 (1991). 29. 0. C. Zienkiewicz, The Finite Element Method, 3rd Edn. McGraw-Hill, New York (I 977).
004s7949/93 $6.00 + 0.00 0 1993 Pergamon Press Ltd
MULTI-STAGE FINITE ELEMENT ALGORITHM FOR EXCAVATION IN ELASTOPLASTIC SOILS E.
COMODROMOS,
TH. HATZIGOGOS and K. PITILAKIS
Department of Civil Engineering, Division of Geotechnical Engineering, Aristotle University of Thessaloniki, Thessaloniki, Greece (Received 12 November
1991)
Abstract-A
multi-stage finite element algorithm for simulating excavation in elastoplastic soils is presented. Based on the conception of the variable domain vector, the non-linear finite element equations are derived in order to take into account the appropriate domain of interest only. Numerical instabilities arising from element and node subtraction during the numerical procedure are overcome by using the static condensation algorithm in the equation solution. Numerical examples demonstrate that the proposed method satisfies the basic criteria of stress-free excavation surface and the principle of uniqueness for linear elastic materials. Furthermore, the notion of the variable domain vector in conjunction with the static condensation algorithm leads to the elimination of any contribution of the removed elements in the numerical procedure and consequently guarantees no violation of the principle of virtual work resulting from the simulation process. In the case of an elastoplastic analysis, a dependency of the provided final solution upon the expanse of the plastified zone has been found. In that case, the one stage analysis is unable to take into account modification arising from topological variations of the structure, such as tangential elastoplastic constitutive matrix. A final solution closer to the exact one is achieved by a multi-stage excavation due to intermediate stages which provides a more accurate behaviour of the previously plastified element.
INTRODUCTION
zero nodal force vector acting on removed nodes, proposed a hybrid method for calculating nodal forces which may preserve uniqueness, but at the expense of added analytical complexity. Borja et al. [13] proposed a variational formulation-based method which in conjunction with the conception of an infinitesimal stiffness satisfies the principle of uniqueness. Finally, Comodromos et al. [14-l 61, have proposed a particular algorithm for nullifying nodal forces arising from excavated elements, which in conjunction with a double pivoting algorithm for the equilibrium solution process satisfies both the principle of uniqueness and the stress-free surface. The simplicity of the conception of the infinitesimal stiffness matrix for the removed elements used by the aforementioned algorithms is offset by the two following disadvantages:
The conception of the early numerical methods for simulating excavations [l-5], was based on the physically coherent and widely adopted hypothesis of a stress-free surface resulting from the governing equations and the boundary conditions of the equilibrium of a continuum body.
An additional criterion of efficiency, accuracy and acceptance known as the principle of uniqueness was postulated by Ishihara [6]. According to the principle of uniqueness, which arises from the principle of superposition, there exists a unique solution independent of the sequence of excavation for a linear, time-independent and elastic material. Ghaboussi and Pecknold [7] demonstrated the validity of the principle by employing a virtual work formulation over the appropriate domain of interest. The theoretically ideal uniqueness is not achieved by many widely used excavation algorithms such as those proposed by Christian and Wong [8], Clough and Mana [9], and Mana [lo, 111. On the contrary, numerical examples of the aforementioned algorithms outlined a significant dependency on the number of excavation stages and discretization of the problem domain. The main conception of the above algorithms consists of limiting the contribution of the excavated elements to the solution procedure by considering an infinitesimal stiffness for the removed elements and trying to nullify the forces acting on removed nodes. Desai and Sargand [12], to overcome the violation of the principle of uniqueness resulting from a non-
1. Violation of the principle of virtual work arising from the contribution of the excavated elements to the equilibrium solution process. 2. The whole procedure leads to a global stiffness matrix presenting infinitesimal pivots corresponding to removed nodes. The use of infinitesimal pivots produces a considerable numerical error which may be overcome by a particular double pivoting algorithm, but at the expense of doubling the solution time. The objective of this paper is to propose an efficient, finite element based, numerical procedure for simulating multi-stage excavation for soils 289
obeying an elastoplastic law with particular focus on the satisfaction against the violation of the principle of virtual work resulting from the simulation of the excavated region. GENERAL
FORMULA~ON
Consider the equilibrium of a body given in Fig, 1, The external forces acting onto the body are surface traction f” and body forces P. The equilibrium solution must satisfy the boundary conditions given by eqns (2) and (3) while at every point within the body eqn (1) must be satisfied too
The presence of a given finite element, e, in the domain of interest, for a given stage, can be arranged by the use of the notion of the variable domain vector (VDV) @,, . The VDV is a vector of scalars which in conjunction with the assemblage process, eqns (5) and (6), determine in a unique way the appropriate domain of interest at a given stage n
~2, = i r=l
Vu(u) -P u=ug
= 0
in f.!
on r,
a(h) = P on f,,
(1) (2) (3)
where V denotes the vector gradient operator, u is the Cauchy stress tensor, up is the prescribed displacement vector, Q is the problem domain, and T,uT, = r are the problem boundaries. In the context of the finite element method, the original problem domain arises from the assemblage of the m finite elements in which the domain has been subdivided. The appropriate domain for a given stage can be arisen from an appropriate assemblage of finite elements. Let e be a finite element of the body given in Fig. 1
12”,= f c= I
[’
dR’ = T@,$, Jn.a
(6)
in which ‘es,, = {fI,&, . . . , S&, . . . ,8&) is transpose of the variational scalar vector with 8:, = 1 when a >n otherwise B:, =0 (Fig. 1) and, Sz, = {Q’, . . . ,ff’, . . . , fy”) is the finite element domain vector. The principle of virtual work is affected by the aforementioned variations and becomes l+z, = YXR>
(7)
where the internal virtual strain energy is given by eqn (g),
The substitution
of eqn (5) to eqn (8) leads to
(4) When applied to excavation problems, both boundaries and domain of interest are stage dependent and eqns (l)-(3) are generalized to allow their variations with excavation stages while eqn (4) is modified in order to take into account the appropriate domain corresponding to a given stage n. The variation of the domain of interest, resulting from the excavation process, involves the subtraction of finite elements according to the stage of the excavation.
Area to be excavated
in which TV’,u’ and @&are the transpose of the strain vector, the stress vector and the variable scalar of element e, respectively. Within the framework of the finite element method, the element strains and displacements for a multi-stage non-linear analysis can be expressed as functions of the global displacements vector ‘II; = N’UI:
(10)
‘6; = B’UZ,
(11)
where N’ denotes the shape function matrix of element e, B’ the stress-strain transformation matrix and Ut the global displacement vector corresponding to stage n and iteration solution k. Furthermore, using the notion of the initial stress method which is convenient for an elastic-plastic analysis 117, 181, the final stresses at kth iteration of stage n, can be expressed as functions of the previous stresses and the actual strains by means of the following equation Fig. I. Schematic variation of domain and boundaries produced by an excavation.
Multi-stage
FE algorithm for excavation in elastoplastic soils which can be rewritten incremental form
in which
ec:=g
the well-known
Kkuk=R” ” ” ex,-R!’ IfI,cr” .
is the consistent tangential matrix obtained by evaluating the variation of the stress vector ‘us with respect to the strain vector ‘.sk If we substitute ind’eqn (9) for the element displacements, strains and stresses, using eqns (IO)-( 12) we obtain
e;,,‘u;
under
(13)
n
Win,= i
291
(18)
The matrix Ki is the global tangential stiffness matrix arising from the assemblage of the tangential element stiffness matrix ki O;,aTBe‘Cy:Be dR’,
(19)
where I” is the residual force vector resulting from the subtraction of the internal force vector R”,, from the external force vector given by eqns (20) and (21), respectively
Been, _ , dR’
e-l s C%
It is obvious that eqn (14) leads to a zero internal virtual strain energy for every element presenting n > a, i.e. when the element is subtracted. The external virtual work for a continuum body with no variations of the domain of interest is given by the well-known equation
+f
0 E,,‘Be ‘a, _ , da’.
R”,,= 2 e=l
TNefSdT’ + F’ (20) e=l I P (21)
Im
NUMERICALIMPLEMENTATION
Applying the notion of the VDV, eqns (5) and (6), and using for the instance F’ and r’ instead of Fj, and r; eqn (15) becomes tl;,, “U; TNefbdR’ P=l s RI
W,, = i
rU;WTNef”dI-e+ ZT;F’.
(16)
It can be seen, that in the case of a subtracted element, the first term of eqn (16) is equal to zero, while the last two are not; that violates the principle of virtual work. However, as will be mentioned later, the component of the external force vector in which the contribution of those terms is stored is eliminated by the static condensation solution process. Finally substituting eqns (14) and (16) to eqn (7) and multiplying from the left by the transpose of the global displacement vector, ZTi, we get the following governing equation
TNef=dT’
(17)
By directly invoking the variable domain vector in the variational formulation the appropriate domain of interest can automatically be obtained without having to renumber nodes and elements of the original mesh. From eqn (19) arises that ‘k, = 0 if a 2 n; the contribution of the initial stress vector became zero in that case as well. However, an equilibrium violation results from the surface forces vector and the concentrated loads of the removed elements. This violation and the mathematical problems arising from the semi-definite form of the global stiffness matrix due to the stiffness matrices of the removed elements can be overcome by the use of the static condensation algorithm in the solution of the equilibrium equations. Static condensation [ 19,201, is employed to reduce the number of element degrees of freedom corresponding to the removed nodes and thus, in effect, to perform part of the solution of the total finite element equilibrium equation system. For each excavation stage an assemblage of the global stiffness matrix and forces vector takes place based on eqn (17) and then a reformation reduction is performed according to the following process: if K,, = 0, then the row and column i of the global stiffness matrix and the row i of the force vector are eliminated as well. The outlined procedure leads to a subdivision of the problem in two independent substructures. The first one concerns the remaining domain and the second one the removed region. The variational equations for the remaining domain are not affected
292
E.
~~ODROM~
by the removed domain and so the problem can be solved with respect to governing equations and modified boundary conditions without any violation. The proposed method has been implemented in a FORTRAN program called nonlinear finite element analysis in geomechanics (NFEAG). It has been verified by comparing predictions with closed-form solutions and other numerical schemes for problems in simulation of excavation sequences. In order to demonstrate the validity of the whole numerical procedure two typical numerical examples are presented. The first, concerns a numerical patch test on a linear elastic one-dimensional excavation problem solved by Clough and Mana [9] and Desai and Sargand 1121. The second one, concerns a twodimensional plane strain excavation problem which has been solved for linear elastic behaviour and an elastic-plastic behaviour with or without hardening. Example 1. Simulation of a linear elastic one-dimensional excavation In order to verify the validity of their algorithm, Clough and Mana [9, 111, performed a linear elastic one-dimensional excavation (Fig. 2) and compared their results with those of the exact solution given below. They removed the upper 20 m by using a oneand eight-element simulation, using eight-node elements. The elastic parameters have been taken as follows: E = lO,O~ t/m2; y = 1.00 t/m2; Y = 0.2 and & = 0.5; the Young’s modulus, the unit weight, the Poisson’s ratio and the at-rest earth pressure, respectively. The same problem was solved by the hybrid method [12], with one-element excavation mesh. The exact solution can be provided by the following linear elastic stress-strain equation
in which f,, and rIij denote the volumet~~ strain and the Kronecker’s delta, respectively.
Area to be excavated
t
A
A
a
Fig. 2. One-dimensional excavation problem.
et al.
Fig. 3. One- and four-stage simulation of the one-dimensional excavation problem.
Taking into account the boundary condition and simulating the excavation by a surface traction of 20 t/m, eqn (22) becomes a
E “=~+(l+v)(l-2v)~~~
vE (23)
from which results a vertical displacement of 3.60001 cm. The same problem was simulated, using the method proposed here, by one and four excavation stages (Fig. 3) by removing one row-element at a time. The solutions arising were identical providing a displacement of 3.60001 cm (exact solution) for the nodes lying on the excavation surface. The vertical displa~ment profiles of the above solutions are given in Fig. 4. It can be seen that the algorithm proposed by Clough and Mana produced a solution highly dependent upon the number of excavation elements. On the contrary, the hybrid method is not affected by the number of excavation elements and stages and provides a solution very close to the exact one. However, both algorithms violate the principles of the one-dimensional form of the problem providing different displacements for nodes lying on the excavation surface. The proposed method, by taking into account the approp~ate domain of interest, eliminates any numerical parasitic error resulting from the cont~bution of the excavated elements and consequently provides the exact solution without any numerical error. Example 2. Simulation of a two-dimensional plane strain excavation The aim of this example is to verify the effeo tiveness of the proposed algorithm in the case of an elastoplastic analysis and furthermore to draw some useful conclusions resulting from such a case. Our interest is focused on an open vertical excavation given in Fig. 5. The mesh is composed of 72 eight-node isoparametric elements and a condition of plane strain is assumed. Nodes along the bottom are fixed in both directions, while nodes along the right-hand and left-hand boundaries may translate freely along the boundaries but are fixed against displacements normal to these
Multi-stage FE algorithm for excavation in elastoplastic soils
293
Drucker-Prager Failure Surface
\< \
/ k
i/ \
/
‘\ \.~
/ ,.:
l. \ ‘.
Fig. 6. Schematic presentation of the cone cap model.
i
Table 1. Summary of parameters used in second example
,,.,‘oneElemcnl . .., . .._+._.A’
Mlnr
1111
Constitutive law Y
Fig. 4. Example 1, vertical displacement excavation surface.
profiles of the
boundaries. In the analyses performed the elements in the upper right-hand corner, bounded by the thick line, were removed in one step and three steps, three row-elements at a time. The aforementioned excavation scheme has been analysed for linear elastic behaviour, perfect elastic-plastic behaviour (Drucker-Prager failure surface), and elastic-plastic behaviour with hardening (cone cap model) (Fig. 6). The parameters used in various analyses, given in Table 1, correspond to the Boston Blue clay [21]. W, D and R are parameters of the elliptical hardening surface, a and k are the Drucker-Prager failure surface parameters and K, G are the bulk and shear moduli, respectively. For elastoplastic analyses the convergence is achieved once the following criteria are satisfied:
Fig. 5. Two-dimensional
@N/m’)
rc,
K (kPa) G (kPa) k (iPa) W D (kPa-‘) R
Linear elastic
DruckerPrager
19.8 0.9 4700 2200 -
19.8 0.9 4700 2200 0.25 10
-
-
CAP 19.8 0.9 4700 2200 0.25 10 0.3 6.1 x 1O-4 4.13
1. Displacement criterion
(24) 2. Residual force criterion
llr2II - II4II < L llrl II ’ f’
(25)
where cdand q/denote the error tolerances which were taken as 10m3 and 10e5, respectively and [I./I is the Euclidean norm.
plane strain open excavation mesh.
294
E.
COMODROMOSet
al.
Linear elastic analysis
Single-stage and three-stage excavation analyses were performed. In the first, the elements 54, 55, 56, 62, 63, 64, 70, 71 and 72 were removed simultaneously while in the second one, a three-step layerby-layer excavation was considered. Figure 7 shows the evolution of the horizontal displacement profile at the vertical boundary of the excavation, line AB on Fig. 5, for the three-stage excavation analysis (configuration R, -+ R, -+ R, + Q,) and the profile corresponding to the single-stage analysis as well (% +RJE). It can be seen that the resulting final solution for both single- and multi-stage analysis are identical, a fact that verifies the principle of uniqueness. Drucker-Prager
= al, + &-k
= 0,
(26)
where 1, denotes the first invariant of stress tensor, J, the second invariant of deviatoric stress tensor while a and k are material constants related to the frictional and cohesive strengths of the material. Strains are elastic for stress paths remaining under the ultimate yield surface, but are elastoplastic if a stress path violates the surface. The same analyses as in linear elastic case were performed; the evolution of the horizontal displacement at the vertical excavation boundary is given in Fig. 8. A small difference between the two final solutions arises from the profile comparison. The
Steo Level
First
Second
Step Level
Third Step Level c
50 mm -
Step Level
Third SteD Level
-
criterion
In this case, an elastic-ideal plastic soil behaviour is considered according to the Drucker-Prager criterion given by the following equation f(a)
First
50
mm
-8-
First
-s
Second
+
Third
Step
+
lnltlsl
PosItIon
-
Step Step
Fig. 8. Example 2, evolution of horizontal induced displacements along the vertical excavation boundary produced by the Drucker-Prager analysis.
maximum difference, observed at node 173, is 1.5% of the displacement corresponding to R, configuration. The plastified integration points at the end of the excavation are given in Fig. 9. Cone cap model
Cap models describe the yielding behaviour of soil with an ultimate yield surface that is fitted with a movable end cap. Both the ultimate yield and cap surfaces are symmetric about the hydrostatic axis. As ultimate failure surface of the cone cap model, the Drucker-Prager surface is used while the strain-hardening elliptical end cap yield surface is described by eqn (27). The movement of the cap is controlled by the hardening function given by eqn (28). Strains are elastic for stress paths remaining within the region bounded by the ultimate yield and cap surface, but are elastic-plastic if a stress path violates any of the two surfaces. The model presented in Fig. 6 consists of a cone shaped Drucker-Prager ultimate surface, fitted with an elliptical end cap, eqns (26) and (27), respectively
./-(I,, ,/&, K) = 111-
L(K)]* + R*J: + [x(K) - L(K)]2
+ +
Flmt
-&
Sacond
(27)
Single-stage K =
Step
& = w[l - eXp(-DX(K)],
(28)
Step
+
Third
Step
+
Inltlal
POIltlon
/
Fig. 7. Example 2, evolution of horizontal induced displacements along the vertical excavation boundary produced by the linear elastic analysis.
where R is the ratio of the major to the minor axis of the elliptic cap, X(K) and L(K) define the intersections of the cap with the I, axis and the failure surface, respectively and K is the hardening parameter, & denotes plastic volumetric strains and W, D are material constants [21,22].
295
Multi-stage FE algorithm for excavation in elastoplastic soils
ii
I I
I
I
I
Fig. 9. Example 2, state of stress at integration points produced by the Drucker-Prager
The horizontal displacement profiles at the vertical excavation boundary corresponding to single- and multi-stage excavation analyses are given in Fig. 10, while the evolution of the plastification zone corresponding to Q,, sl, and i& con~gurations are given in Figs 1l-13. The difference between horizontal displacement at node 173 provided by single- and multi-stage excavation is 2.5% of the total displacement. Remarks
A comparison between the aforementioned displacement profiles for linear elastic analysis, the Drucker-Prager elastic-ideal plastic analyses and the cone cap model is presented in Pig. 14. Figures 15 and 16 show the displacement history of nodes 168, 173 obtained from the single-stage excavation (Q,, --) C&) multi-stage and the three-step excavation (Q, -+ a, -+ a2 -+ 12,). It can be seen that the cone cap model gives the greater displacement and the greater difference between single- and multi-stage analyses. This is due to the expense of the plastified zone arising from the delimitation of the linear elastic zone by the elliptical end cap. The accuracy and the capability of the proposed method for solving non-linear excavation problems with particular application to elastoplastic soils with or without hardening outline its effectiveness and usefulness.
instabilities arising from element and node subtraction during numerical procedure are overcome by the use of the static condensation algorithm in the equation solution. Numerical examples demonstrated that the proposed method satisfies the basic criteria of stress-free excavation surface and the principle of uniqueness. Furthermore, the notion of variable domain vector in conjunction with the static condensation algorithm eliminates any contribution of the removed eiements in the numerical procedure and consequently guarantees no violation of the principle of virtual work resulting from the simulation process.
A multi-stage finite element algorithm for simulating excavation in elastoplastic soils has been presented. Based on the conception of the variable domain vector, the non-linear finite element equations are derived in order to take into account only the appropriate domain of interest. Numerical
Step Level
First
bird Step Level -
+ CONCLUSIONS
analysis.
50 mm -
Singlestage Step
-8-
Flrllt
+-
Second
-g
Third
Step
-4%
lnltlsl
POIItlon
Step
Fig. 10. Example 2, evolution of horizontal induced displacements along the vertical excavation boundary produced by the cone cap analysis.
gelastic mc,
~Cotncar
~Tenslm
[‘ZiaFailure
Fig. 11. Example 2, state of stress at integration points produced by the cone cap analysis at stage one.
Fig. 12. Example 2, state of stress at integration points produced by the cone cap analysis at stage two.
Fig. 13. Example 2, state of stress at integration points produced by the cone cap analysis at stage three. 296
Multi-stage FE algorithm for excavation in elastoplastic soils
297
modification arising from topological variations of the structure, such as tangential elastoplastic consti-
tutive matrix. A final solution closer to the exact one is achieved by a multi-stage excavation due to intermediate stages which provides a more accurate behaviour of the previously plastified element.
50
-yt
LItmar-Elaatlc
+-
Druckw-Prager
-a-
conm Cap
-+
lnltlsl
REFERENCES
mm
Pomltlon
Fig. 14. Example 2, comparison between the linear elastic, the Drucker-Prager and the cone cap final solutions.
0
1
L
3
4
Horizontal
5
6
Dfsplacement
7
6
(cm)
Fig. 15. Comparison between the displacements of the node 168 produced by the linear elastic, the Drucker-Prager and the cone cap analysis.
-1
4
0
1
tioYizon:ai D:splafxmk7t
ikn)
0
9
IO
Fig. 16. Comparison between displacements of the node 173 produced by the linear elastic, the Drucker-Prager and the cone cap analysis.
In the case of an elastoplastic analysis, a dependency of the provided final solution upon the expanse of the plastified zone has been found. In that case, the one-stage analysis is unable to take into account CAS 46124
1. J. M. Duncan and C. Y. Chang, Nonlinear analysis of stress and strain in soils. J. Soil Mech. Found. Div., ASCE 96, 1629-1653 (1970). 2. J. M. Duncan and P. Dunlop, Slopes in stiff fissured clay and shales. J. Soil Mech. Found. Div., ASCE 95, 81-104 (1969). 3. G. W. Clough and J. M. Duncan, Finite element analyses of Port Allen and Old River locks. University of California, Berkeley, Report TE-69-3 (1969). 4. G. W. Clough and J. M. Duncan, Finite element analyses of retaining wall behavior. J. Soil Mech. Found. Div., ASCE 97, No. SM12 (1971). 5. P. Dunlop and J. M. Duncan, Development of failure around excavated slopes. J. Soil Mech. Found. Div., ASCE 96, 471-493 (1970). 6. K. lshihara, Relations between process of cutting and uniqueness of solutions. Soils Foundations 10, SO-65 (1970). 7. J. Ghaboussi and D. A. Pecknold, Incremental finite element analysis of geometrically altered structures. Int. J. Numer. Meth. Engng 20, 2051-2064 (1984). 8. J. T. Christian and I. H. Wong, Errors in simulation of excavation in elastic media by finite elements. Soil Foundations 13, l-10 (1973). 9. G. W. Clough and A. I. Mana, Lessons learned in finite element analysis of temporary excavations. In Proc. Zhd Int. Conf on Numer. Meth. in Geomech., ASCE, Blacksburg, VA (1976). 10. A. I. Mana and C. W. Clough, Prediction of movements for braced cuts in clay. J. Geotech. Engng Div., ASCE 107, 759-777 (1981). 11. A. I. Mana, Finite element analysis of deep excavation behavior in soft clay. Ph.D. thesis, Stanford University, Stanford, CA (1978). 12. C. S. Desai and S. Sargand, Hydrid FE procedure for soil-structure interaction. J. Geotech. Engng, ASCE 110, 473-486 (1984). 13. R. I. Borja, S. R. Lee and R. B. Seed, Numerical simulation of excavation in elastoplastic soils. Inf. J. Numer. Anal. Meth. in Geomech. 13, 231-249 (1989). 14. E. Comodromos, T. Hatzigogos and K. Pitilakis, Finite element excavation in elastoplastic soils. In Proc. &d Eur. Spec. Conf: Numer. Meth. Geotech. Engng, pp. 573-584. Santander (1990). 15. E. Comodromos, A finite element algorithm for simulating excavations in elastoplastic soils. 5rh Young Geotech. Eng. Co& Grenoble (1991). 16. E. Comodromos, Contribution to the analysis of excavation by the finite element method. Ph.D. thesis, Aristotle University of Thessaloniki (1991). 17. 0. C. Zienkiewicz, S. Valliapan and I. P. King, Elastoplastic solutions of engineering problems-initial stress finite element method. Int. J. Numer. Meth. Engng 1, No. 1 (1969). 18. C. S. Desai and F. J. Abel, Introduction to the Finite Element Method. A Numerical Method for Engineering Analysis. Van Nostrand Reinhold, New York (1972). 19. K. J. Bathe and E. L. Wilson, Numerical Methoak in Finite Element Analysis. Prentice-Hall, Englewood Cliffs, NJ (1976). 20. E. L. Wilson, The static condensation algorithm. Int. J. Numer. Meth. Engng 8, 199-203 (1974).
298
E. CO~OD~O~~ et al.
21. D. N. Humphrey, Design of reinforced embankments.
Ph.D. Thesis. School of Civil Engineerinn, Perdue University, West Lafayette, IN (1986). -. 22. W. F. Chen and G. Y. Baladi, Soil Plasticity-Theory and Implementation. Elsevier, New York (1986). 23. C. Y. Chang and J. M. Duncan, Analysis of soil movement around a deep excavation. J. Soil Mech. Found. Diu., ASCE 96, 1655-1681 (1970). 24. V. S. Chandrasekaran and J. W. King, Simulation of excavation using finite elements. Tech. Note, J. Soil Mech. Found. Div., AXE 100, 1086-1089 (1974).
25. C. S. Desai, Soil-structure interaction and simulation problems. In Finite Element in Geomechanics (Edited by G. Gudehus), pp. 209-250. John Wiley (1977).
26. C. S. Desai and H. J. Siriwardane, A concept of correction functions to account for non-associative characteristics of geologic media. Int. J. Numer. Anal. Meth. Geomech. 4, 377-387 (1980).
27. G. Dhatt and G. Touzot, In Une Prbentation de la MPthode des Uments Finis, Znd Edn (Edited by S. A. Maloine). Collection de I’Univercit~ de Compiegne, Paris (1984). 28. D. F. Stolle, An interpretation of initial stress and strain methods, and numerical stability. Int. J. Numer. Meth. Engng 15,399-416 (1991). 29. 0. C. Zienkiewicz, The Finite Element Method, 3rd Edn. McGraw-Hill, New York (I 977).