- Email: [email protected]
Interaction problems in soil-structure mechanics with material nonlinearity
INTERACTION PROBLEMS IN SOIL-STRUCTURE MECHANICS WITH MATERIAL NONLINEARITY KENNETHRUNESSON,HARALDT~~GNFORS and NILS-ERIK WIBERG Chalmers University of Technology, Department of Structural Mechanics, Giiteborg, Sweden (Received 2 October 1979;received for publication 29 January 1980) Abstract-Finite element analysis of soil-structure stress interaction and sand-clay flow interaction is discussed in the paper. The soil may be a one- or two-phase medium, depending on if fully drained, partly drained
(consolidation)or undrained conditionis considered.Transient behaviouris investigated. Yield hinge theory is adopted for the structural members while visco-plasticity, plasticity or a combination of visc&plasticity and plasticity are adopted for the soil skeleton to account for inviscid yield and creep. Adopted yield criterions are of cone type and according to the Critical State Theory. Finite element calculations for the stability analysis of a natural slope in clay and the interaction of a frame resting on a consolidating elastic-plastic subsoil are presented. 1. INTRODUCTION
Many problems include subregions with quite different physical behaviour. However, they interact and give mutual dependence, and it is of utmost importance to account for this in a unified analysis. Soil-structure interaction on structural level is described for instance by Goschy[l], and for a single structural member by Davidsson and Chen[2]. In Zienkiewicz et a/.[31 fluidstructure dynamic interaction is described. In interaction problems, substructuring technique is very convenient both for the problem description and for the numerical calculation. It is convenient that each substructure contains a special problem based on a physical law. The special problems (substructures) are coupled by interfaces by which special physical phenomena can be handled, such as crack, slide, pore pressure discontinuity etc. In this paper the interaction problems in soil due to transient flow are investigated. The paper starts with a review of some general theory for saturated soil which has to be treated as a two-phase material. The motion of the two phases (the soil skeleton and the pore water) is coupled due to interactive forces. The basic mechanical relations for the coupled problem are given elsewhere, e.g. Sandhu[4], Runesson[S]. In judging the magnitude and progress of the settlements of a soil layer, the choice of a reliable constitutive model for the soil skeleton is of outstanding importance. For soft clay it seems appropriate to include plasticity in the constitutive equations. Creep effects are important in many Swedish clays, especially when the organic contents are significant. There seems to be experimental evidence for the use of visco-plasticity as observed from undrained triaxial tests, Larsson[6]. However, though
the creep rate has been tested under fairly restrictive stress conditions, use will be made in this paper of an associated creep law according to the theory by PerzynaV]. One type of interaction is the interaction between layers of silt (or sand) and the embedding clay. Because of the larger permeability in the silt (sand) layer, the pore pressure may increase, for example due to heavy raining, much more rapidly in the silt than in the clay. As a result the effective stresses decrease first in the silt and successively in the clay. thus reducing the strength. A number of landslides in natural slopes in Sweden during recent years have made this problem paramount, Fig. la. Interaction of different types of soils claims for the use of different types of constitutive properties (e.g. yield conditions) in the same problem. In most construction problems in practice it is not sufficient to analyse the superstructure and the soil separately. In general it is more important to account for true interaction in the case of a slab-wall structure than a beam-column structure, which is more flexible. The redistribution of internal forces caused by the interaction phenomena with the soil is time dependent (consolidation and creep). Long term effects are extremely important in gravity structures. In Sweden many aquifers consist of fairly thin sand layers, which are confined from above by a thick deposit of soft clay and from below by solid rock. An important consequence of a disturbance of the hydraulic equilibrium is the occurrence of pronounced settlements at the ground surface. It is obvious that a complete analysis must include both the transient see page in the aquifer itself as well as the more delayed settlements and pore
’ settlements
clay
rock
Fig. 1. Interaction problems. 581
^ ” \ tunnel
582
KENNETHRUNESSON et al.
water flow in the clay and, finally, the behaviour of the structure under the time dependent settlements. Figure lb shows a typical situation where a frame interacting with the soil is subjected to hydrodynamic loading from leakage into a rock tunnel in the bedrock below. 2. MECHANICS FOR SATURATED SOIL
2.1 Effective stress hypothesis Saturated soil may be idealized by a mechanical model consisting of basically two phases, a solid phase (the soil skeleton) and a fluid phase (the pore water), mixed together. According to the effective stress principle by Terzaghi[8] the total stress &ii is equal to the difference between the effective stress ait in the soil skeleton and the pore water pressure p (tension stresses positive). Crij = Uij - p&j*
(1)
It is to be noticed that stresses are measured per bulk area. By definition all deformations in the soil are coupled to the effective stresses.
2.2 field equations Adopting Cartesian c~rdinates we may summarize the basic relations for saturated soil. Before any further approximations are introduced, the equations are valid both for sand and clay. ~uiiib~um equations: - (Fij,j - Uij)
t p,i = 0
(2)
meability tensor, and h,, are functions expressing the viscoplastic strain rate. Explicit time hardening has to be contained in hki for an adequate description of the behaviour of many Swedish clays. 2.3 Boundary conditions The boundary of a finite body of saturated soil is divided into two parts in two different ways. On disjunct parts relevant boundary conditions have to be imposed on displacements ui and tractions Ti as well as on the excess pore pressure p and the draining flow Q.
3. NON-LINEAR CONSTITUTION OF THE SOIL SKELETON
3.1 General Because of the great variety in ~haviour of soil depending on sedimentation environments, porosity, water contents, organic contents, cementing interactive forces, etc., there is little possibility to unify the behaviour of soil in one single constitutive model. In particular, many Swedish clays are soft with high porosity and contractant under deviatoric loading beyond the preconsolidation stresses. Even under rapid, undrained loading a dramatic yield is obtained at the state of preconsolidation, S~liforsr9]. In some cases the behaviour is almost perfectly plastic. In addition, creep effects are significant in the overconsolidated region, even for moderate additional stresses. These properties indicate that the choice of an elastic-viscop~asti~-plastic model is well motivated. Pure creep and inviscid plasticity are obtained as special cases of the more general model.
where
The weight of the skeleton represented by the initial volume forces II? is equilibrated by the initial stresses &. Additions loads are applied only on the boundary, thus making (2) homogeneous. Further, p is the excess pore pressure (diffusive pressure), since the hydrostatic pressure equilibrates the pore water bulk weight. Continuity condition: & + 4k.k = 0 (3) where rkk is the volumetric strain rate and qi are the pore water velocities relative to the soil skeleton. A dot denotes differentiation with respect to time. Equation (3) arises from the assumption of the pore water being completely incompressible. Kinematic equations (small strains): Cij= !$U,.j + Uj.i), qi = - P,i
(4)
where displacements are denoted by tci,strains by +, and diffusive forces by gi. Constitutive equations for the soil skeleton tincluding elastic-viscoplastic-plastic behaviour) and for the pore water flow (Darcy’s law): Cij = S&(&f - hki), (3 I s$, = S$&%z,,, &), hkj = h&&n, tf qi = Kigi
(6)
where S& are components of the efasti~-plastic tangential stiffness tensor, Kii are components of the per-
3.2 Elastic-viscoplastic-piastic modes The rheologicat model for the material model used is given in Fig. 2a. Associated with the two sliders are two families of closed surfaces in stress space, Fig. 2b, the inviscid yield surface (a = I) and the quasistatic yield surface (a = 2), defined by f~=O,f~=f~(c+ii,K,),(y=1,2
(7)
where K* are h~deni~ (softeni~) parameters. In the present application the inviscid yield surface f, = 0 is assumed to be hardening (softening). The plastic part of the volumetric strain is chosen as the hardening measure K$,i.e. restriction is made to pure density hardening. The quasistatic yield surface f2 = 0 is assumed to be perfectly plastic so that ~~ is irrelevant. The functions fa are similarly defined such that the surfaces f_ = 0, i = 1,2, completely coincide for a certain amount of softening in f, = 0. While f, I 0 must always hold, there is no such restriction on &. Thus, f:! = c, c > 0, define subsequent loading surfaces. Depending on the position of the stress point in stress space different characteristic stress states are obtained, Fig. 2b. The kinematic decomposition of the strain rate is (small strains) g=~;+&+&
(8)
Using associated flow rules (which may be questionable for certain choices of fa) both for the plastic part, i&, and for the creep part, i;, in accordance with Perzyna[7], we obtain the consti~tive law in Stiffness form according to (5) with the explicit expressions for
583
Interaction problems in soil structure mechanicswith materialnonlinearity 412
T
, elastic-wscoplastlc-plastic elastic-wscoplastlc
k, 2 k2 (4
(b) Fig. 2. (a) ID-rheological model, (b) loading surfaces in stress space.
the tangential stiffness (9) where
D=af,sf, afl_
-_afl afl &Yij "kra~krawk) aukk
and for the creep rate
hij
=
/I I
v(t)iO(F)>$g 1,
3.4 Hardening Based on experimental observations the material behaviour of silt and sand can be described as shear strain hardening of the extended von Mises yield criterion according to (11) with no cohesion, i.e. pb = 0. Initial yield corresponds to a cone with the slope M”’ = M., Fig. 4a. The material hardens in a drained biaxial test as in Fig. 4b, and the ultimate stage is defined as a cone with the slope M,. Intermediate cones are defined by the slope M as a function of the deviatoric plastic strain l Qp. (13)
(10)
8,
where y and @ are creep functions, F is a non-dimensional variable obtained by scaling f2 with a reference parameter. The length of a second order tensor ‘I7 is denoted ]&I = (Tk,Tk,)“*. Further, for a SC&r X, (X)=X forxrO,(x)=Oforx
(11)
containing extended von Mises (compression cone, extension cone), Drucker-Prager, etc., Fig. 3a. Surfaces associated with the Critical State (CS) theory are defined by, Fig. 3b, f~=4*tM*@,-p6)@)-p~),a=1,2
where fi(e/)=
-&)i)
(14)
with the notations ti(eqP) = 4(eq%rb M M’=&j,M,=
-
M, z- 3GeYp 1- M,/3’ 2~0 ’
The relation between the deviatoric plastic strain E,” and the deviatoric stress 4 in a drained biaxial test is shown in Fig. 4b, where the non-linear part is assumed to be a hyperbola. The test starts at the isotropic consolidation stress pb. The shear modulus in the elastic region is G. A hardening rule for clay suggested by Janbu[l I] is used in combination with yield surfaces associated with the Critical State theory, (12),
(12)
where, in particular, p)) is the isotropic preconsolidation pressure, pb is a cohesion parameter, and M is the slope of the Critical State Line (CSL), defined as h = q-M(p’-p6)=0. An anisotropic extension of (12) has been developed recently, Runesson and Axelsson[ 101,Runesson[S].
Ii& --Et X&*+2(&
PC
‘(I)_
-Pro
’ (1)
(15)
ev(kmqP)
where p :j is the initial preconsolidation pressure and k,,, is a nondimensional constant. The choice seems reasonable for Swedish soft clays. 3.5 Creep Recent tests (undrained) on Swedish clays, Larsson[6],
Q
,h=O /
,&
I
kSU , anisotropic
Fig. 3. Isotropic yield criterions: (a) Extended von Mises etc., (b) Critical State.
locus
KENNETH RUNESSON et al.
584
(a)
(b) Fig. 4. Yield criterion for silt. indicate a linear creep law in the overconsolidated region
while an exponential dependence similar to the one by Singh and Mitchell[l21 seems reliable for stresses in the post-yield region. In addition, we have to take explicit time hardening into account. Details about this can be found in Runesson[S]. 4. UNIFIED NUMERICAL ANALYSIS
4. I General A unified finite element analysis is devised for the soil-structure interaction problems discussed above. Restriction is made to ZD-cases. One very important problem is the stability of natural slopes, Fig. 5. In Sweden these problems are of current interest because of recent landslides in high sensitive plastic clays, Maripuu and Berntsson[l3], have caused very severe consequences like destroyed buildings and lost lives. The problems have a 3D-character, but in many cases it seems reasonable from the physical point of view, and necessary from the computational point of view, to approximate the problem as a 2D plane strain problem. Plasticity and creep together with the magnitude of the pore water pressure play an important role for the overall behaviour. Because of the softening character of the clay, local collapse of the clay may cause a progressive collapse. The local collapse may be caused, for example, by creep or increase of the pore water pressure. The small “local superstructure” has very little influence on the total performance of the natural slope, Fig. Sb.
It is important to analyze the influence of rapid change of the pore water pressure in the silt or sand layers embedded in the clay, Fig. 5b. Increase of pressure may be caused by a period of heavy raining. In the example in Section 6.1 a preliminary calculation based on a one phase material model is shown. Another important problem concerns the long term, time-dependent behaviour of thick clay deposits due to consolidation and creep, which cause settlements and redistribution of stresses in interacting structures. As a result severe damages on buildings may occur. Figure lb shows a situation where there is a groundwater leakage into the rock tunnel from the sand layer. The resulting pore pressure drop in the sand layer will, however, very soon become stationary. This decrease of pore pressure in the sand layer starts a consolidation process in the clay layer, which results in settlements at the ground surface. The frame interacts with the soil so that the settlements cause additional displacements and stresses in the frame. A third case is to analyse the swelling problem. In many cities a lowering of the groundwater table has led to consolidation and thus settlements. By infiltration it is possible not only to stop the consolidation but also induce the surface to raise again to a certain extent. The calculation is valid only for a moderate increase in the pore pressure. The problem is nonlinear and time dependent. The universal parameter is the time t. Loads W, boundary conditions and constitutive parameters may be time dependent, Fig. 6a.
i Lic+Ls
Fig. 5. (a) 3D-natural slope (valley), (b) sdt layer in natural slope, ZD-plane strain.
W
w,
‘,
_____--
I
(a)
\ W,i/i, Fig. 6. Real- and fictitious time.
(b)
585
Interaction problemsin soil structure mechanics with material nonlinearity In order to handle an instantaneous load a fictitious time parameter i is defined. The instantaneous load is applied to the structure as W = W,di, at the real time t,, Fig. 6b. Thus, the action spectra can be applied to the structure in the time space (t, t), see Fig. 6c. 4.2 FE-discretization in space
A discretization in space of an appropriate variational formulation leads to the semi-discrete Galerkin method for problems with time dependence. A system of ordinary differential equations in time is obtained. The space discretization used in the computer program GEOFEM [ 141,for calculation of the 2D-problems shown below, is based on a parametric finite element with a biquadratic displacement approximation and a bilinear excess pore pressure approximation in local coordinates. The element geometry is mapped from quadratic regions in the local coordinates by biquadratic functions. For coupled consolidation such an element may be termed quasi-isoparametric, while for the pseudo-theory (ground water flow analysis) the element will be superparametric.
terms have the required regularity to match the theoretical order of accuracy. As the coefficient matrix will be constant during a time-interval, an extra stage in the method corresponding to the estimation of the truncation error only requires an extra back-substitution provided an elimination solution method is used, which is the case in the computer program GEOFEM. The particular method adopted is called (2,4)-W. It is a 2-stage method of order 2. Applied to the problem
where f is a column matrix of the same dimension as j, the method is defined by j,=&&3k*) W(& 0, P)g, =f(&> to) W(At, (Y,P)k,=f($,t;AtE,,
4.3 Numerical time integration The non-linear initial value problem arising from the FE-discretization of the consolidation problem can be formulated in matrix form as _d6 A@&=-B6t
IC(O)=co.
V(t)tA(& t)
WI
The nodal values of the displacement S and the excess pores pressure p’ are collected in the matrix 5, while the applied load U diffecentiated with respect to time and the known drainage P are collected in V. A pseudoload due to creep, f?‘(tj, t), is contained in A.
The square matrix A, which is state dependent, and B, which is constant, are defined by
t,,t;At
(18) -;AfaPk,
W(At, (Y,P)=Z-AtaP where P is an arbitrary real quadratic matrix and a is the inverse of the largest root to the second Laguerre-polynomial, i.e. a = 1 - l/6(2). The method has a built-in local error estimate of third order. Details may be found in[15]. Due to the certain kind of non-linearity introduced in (16) at the adoption of an elastic-viscoplastic-plastic model, some modification is made at the practical application to the problem of consolidation. This is motivated from a solution economy point of view. The main modification is that the tangent stiffness matrix S is evaluated for t& only and is kept unchanged through all stages during a single time step. Analogeously to solid analysis, the fulfilment of the yield condition in the end of the time interval results in a lack of equilibrium, which is compensated for by applying a fictitious residual load. This load gives rise to additional displacements and pressures. An iterative procedure leads to the simultaneous satisfaction of the yield condition and the variational equilibrium equation. 5COMPUTER
PROGRAM GEOFEM
General A computer program GEOFEM [ 141for the solution of non-linear and time-dependent geomechanical problems is under development. Main characteristics are: 0 Simplicity and generality in the program structure 0 New types of problems can be added in a simple way 0 Input data is given free format introduced by key words 0 The calculations are steered by simple input commands 0 Substructuring technique is used 5.1
where S is the tangential stiffness matrix for the skeleton and C is a rectangular coupling matrix. Further, K is the permeability matrix and GS is associated with a convective flow boundary condition. In the program GEOFEM a single-step method of linear semiexplicit Runge-Kutta type is used to trace the solution in time. The method used is basically a so called W-method, a class of methods which has recently been proposed by Wolfbrandt [151. W-methods are convenient for stiff systems like (16). An automatic prediction of the next time step is based on the calculation of the principal error term, which is determined by a comparison of the results obtained from two methods of different order (imbedding technique). Hence, for small time steps a reliable estimate for the local truncation error is possible, provided the non-linear
5.2 Substructuring. Interfaces The structure is described by a number of substructures coupled by interfaces. The interfaces may either be a logical coupling, or a coupling containing stiffness (or other material) properties. By use of the interfaces cracks, friction, sliding, discontinuities in pore pressure,
KENNETH RUNPSONel 01.
586
and so on can be described. Each substructure has a simple description and the output is obtained separately for each substructure. The substructures may have different mesh sizes, which means that master nodes are used at the interfaces in order to fit the different meshes together. Two cases of interaction from soil-structure mechanics are shown in Fig. 7a and b. Master variables li” in the master nodes are introduced by the transformation (ti refers to a certain substructure)
The normal stresses in situ in the xyz-system (~~~,u~~,u~~) are assumed to be the principal stresses. The total vertical stresses and pore pressures are calculated as follows, Fig. 9:
a=o= i p. =
g=fi*“*
(19)
The stiffness relation for interface nodes are obtained by changing the stiffness relation 5% = U to D’SDC” = o”, Ljm=Dto.
(20)
These multiplications are made in condensed form. Addition of stiffness properties in the interfaces yield the final equations
[
sl;,Ys2;s;][ ;;:I=[;;:]
(21)
where index 1 and 2 refers to substructure 1 and 2, and I refers to the interface. 6.~CAL 6.1.Stability analysis
EXULT
6.1.1 Problem. Basic data. A stability analysis of a natural slope from Tuve in GGteborg, Sweden[l3], is performed. The stope is shown in Fig. 8. We want to estimate the safety factor for the slope. The safety factor is here defined as the factor fi, which multiplied with the total weight of the slope gives instability. The influence of a silt layer on the safety factor is studied.
2 E LZ&zrl z E [z,, &I z E [Z%z,l
P&z -Lo), PO& - zo)+ PC?& - z,), PO&, - zo)+Plzg(&-~I),
pgtz - z2), .?E [zo, ZJ apug@ - z,) + APO,p,,g = 10kNlm’, z E IZI, zzl z E I&, ZJ. i %vg(z - ZI),
The pore pressure p. is dependent of the flow conditions, because of an inclining ground water level. This is accounted for by the factor a. For pure hydrostatic conditions a = 1. Figure 9a shows the distribution of the total vertical stress aVo and the pore pressure p. at different depths, where Ape means a drastic increase of pore pressure in the silt layer (which may be caused by season variations). The local increase of the pore pressure in the silt layer may be approximately calculated by simple transient aquifer analysis. This flow problem becomes stationary in a short time (some weeks or a couple of months). The failure will not take place in the silt layer but in the clay nearby, which justifies the calculations made here. The effective stresses in situ will thus be cr:o=a;o-PO,
(22)
(aLo= -aLo)
0:o = KoaLo, (&o= -Cr:o= -&o). 6.1.2 Used method. The analysis is here performed as an undrained calculation for a one phase material (total stress analysis). For an isotropic elastic material the effective mean pressure p’ is then constant during the / frame
f ting
sltt layer master node (Interface with physical
la)
(b)
Fig. 7. Interfaceswith master nodes.
r
+...__._
-_,x
i z
12m
----_.-._-..~-.____
250m
__
I
Fig. 8. Natural slope
in
Tuve.
4
Interaction
problems in soil-structure mechanics with material nonlinearity
587 kN/m
50
I,(X)
z,(x)
I
1
”
”
,,
@‘“=I5
1
t
9.2
i
4
P0,g=17.5kN/m2 I
z(m)
04
(4
Fig. 9. Distribution of: (a) vertical stresses Q, and pore pressure p0 and, (b) undrained cohesion c, for a vertical profile. loading. The strength can be represented with a deviatoric parameter which is a function of the effective mean stress in situ. Von Mises yield condition is chosen which is a special case of the extended von Mises (DruckerPrager) with
responds to the true deviator strength q. = 2c, for the initial stress state, is determined from (25) by setting @, = 0, Fig. IO. The condition qd = q. for p’ = PAgives, if sin & < 1, c,, = pb sin & t cd cos &,
c$,,= 0.
(26)
C#J = r#J”= 0, c = c,. The determination of the undrained strength c, is demonstrated below. For synchronic loading the proportionality factor serves as a fictitious time variable. The nonlinear initial value problem from the FE-descritization can be formulated in matrix form as (a part of (16))
The effective mean stress pb in situ is defined by, the state stress (TX0 = uyo = (aLo, o:o, uL0) with K do, ( do = - dd as p;=
_1/ it/c = 3
;(I+2K&r:o.
(27)
For clay it is also possible to calculate the C.-value from the undrained shear strength 7fU obtained from (23) vane borer tests. with the constraints of weak equilibrium
&(a - u,,) = Uf, b, =
I ”
cu =
A&‘o d V
(24)
and the yield condition. 6.1.3 Calculation of c,. The average cone yield condition (concerning effective stresses) can be formulated as 18sin&, , 18co&, (25) qd = 9_sinz,$dp ‘m ‘d. The undrained
cohesion parameter
c., which cor-
;[(1-
Ko)*a:a + 3&"*.
It is then assumed that the normal stresses in situ remain unchanged during the test. For the problem in Fig. 8 with the ground water table at the ground surface, the data according to Fig. lob have been chosen. Inserted in (26) it gives the undrained cohesion c. shown in Fig. 9b. A second calculation is made for an increase Appoof the pore pressure, which gives the change in the c.-value AC. = f( 1t 2K,,)( - Ap,,)sin c#+, 6.1.4 Numerical
calculations.
Figure lla shows the
4
T
ApI
Clay:C,,=IOkN/m’ Ko= 06 &= 20” a = I.1
Sand:
(28)
cd=0 K,J=0.6 &=36’ Y 50.495 (b)
P’
(a)
Fig. IO. (a) Simulation of yield condition with undrained parameters, (b) soil data.
KENNETH RUNESSON ef al.
/?=I.55
1.75
1.55
2.15
1.95
2 15
(b)
(a) Fig. 11. Slope: (a) displacement picture for p = 1.65(displacements magnified I5 times), (b) gowth of plastic zone.
finite element mesh and the displacement picture for the load factor p = 1.65. Instability is obtained for the load factor 2.15. The growth of the plastic zone is shown in Fig. lib. Figure 12 shows the displacement picture and growth of the plastic zone for the case with a drastic pore pressure increase of 40 kN/m* in the silt layer giving AC. = - 10kN/m*. Instability in this case is found for the load factor 1.7. In particular, we observe the large shear deformations in the silt layer and the development of the “circular” yield zone. A local zone shown in Fig. 12 (as a shadowed area) with a low shear strength value diminishes the load factor to 1.4. We observe that the failure load is very sensitive to local areas with low shear strength. A more sophisticated and true analysis ought to include creep and softening behaviour of the soil.
6.2 Soil-frame stress interaction and sand-clay fZow interaction The leakage problem in Fig. 13 is studied under consideration of different material properties and different loadings. Both the flow interaction between the sand and the clay layer and the stress interaction between the soil and the frame have to be accounted for. The long term effects from consolidation and creep are studied. The soil displacements and the displacements and redistributed forces in the frame are obtained as results. 6.2.1 Basic data. The soil profile in Fig. 14a is considered under two different loading alternatives. The initial pore pressure distribution is assumed to be hydrostatic. The vertical effective pressure (compression stress) in
Fig. 12. (a) Displacement picture for p= 1.65(displacement magnified 15 times), (b) growth of the plastic zone.
589
Interaction problems in soil-structure mechanics with material nonlinearity
sym
1
7m
*
frame,
c
Sand
o .
- k=lO-s
.
II
I
7m
I
T
I I
I I
.,* II
.
~
III.
o
. .
.
3m
Ah,hhhh
IOm
3m
3m
Fig. 13. Clay-sand flow interaction and frame-soil interaction.
IOm Soft
clay
Fig. 14. (a) Soil profile, (b) initial stress a:,, and preconsolidation stress 0: distribution with depth in situ is assumed as a:o=&,
p’ = Pm - PW.
A normal value for Swedish clay is p,,,g=16kN/m3 giving p’g =6 kN/m’. It is assumed that the preconsolidation pressure ail), which is the equivalent to ~21)for the anisotropic extension of (12),[5], can be estimated as -‘(I) = a:( 1+ 2&)/3 PC where, in the present example, Z&=0.8 is chosen. The distribution of a:, which is the vertical preconsolidation stress obtained at an oedometer test, is assumed as, Fig. 14b, I 2.0~~kN/m*, 01z12m (Tc=I za kN/m*, 2~z~lOm where cx is a constant. The constant value of (T: in the upper part of the clay layer corresponds to a dry crust. The quasi-static preconsolidation pressure & is assumed to be equal to p& the mean effective pressure in situ. It can also be calculated as B:“‘/OCR
where the in situ value of the overconsolidation ratio OCR is OCR+= I
2.OiYlp’gz, 0 5 z I2m ( fflp’g, 2 5 z I 1Om
Beneath the dry crust a slight degree of overconsolidation is assumed by the choice of cr=6.6 (OCR= 1.1). The elastic constants used are E=3000 kN/m’ and v=O.30. The internal friction angle is a=30 (giving M) and a small cohesion is assumed, c=S kN/m* (giving pb). To define the hardening along the normal consolidation line at &-consolidation, the value k,,,=4 has been chosen. The isotropic permeabiiities for the clay, sand, rock and cracked rock are found in Fig. 13. 6.2.2 Numerical calculations. The finite element mesh is shown in Fig. 15. The figure also shows the different types of elements; beam elements (3 variables at each node), consolidation elements (displacements li and pore pressure @), and ground water flow elements (pore pressure 6). The drilling of a rock tunnel initiates a ground water leakage from the sand layer through the zone of cracked
KENNETH RUNESSON et al.
SnuoP%sY”re
I
, soll-structure stress interaction
4 ’
r-_A clay layer
I id -~- ~1
slow transient flow
11
, sand-clay flow interactlon
! P--i w 4 sand and rock, rapid + c transient flow
I I
j
I
1 +I3
___ It-L Fig. 15. Finite element mesh. s
n
clay
t--o
. .. sand ,. . hAx
pore pressure proflle stationary solutlon for the decrease ,of the pore pressure at the upper vurface of the sand layer pore pressure 1 ,
rock symmetry 0x15 after 0.02years
\
hydrostatic
\ \J
/
-
\
: _t,t_+_. 100
3
\
:
:
I
pressure
:
:
,pkNh2
4
300
200
100
tunnel
p kN/m’
Fig. 16. Pore pressure protiles at t = 0.02yr (no creep, no surface load). W
,t tunnel drilling
’ mainly flow in sand layer
consolidation external loadmg
Fig. 17. Loading procedure in real and tictitious time.
591
Interaction problems in soil-structure mechanics with material nonlinearity
displacements O.lJm
Fig. 18. Displacements at (a) t = 0.5 yr and (b) t = 10yr (no creep, no surface load). t=O.Oiyr
sym
I
Fii. 19. Spread of yield zone at t = 0.5 yr and t = 10yr (no creep, no surface load). Section 20
40
Section II -11
I-I
60
60
too :
, p(kN/mz)
\.
k
,=aa
6
Depth (ml
t = IOyr
=IOyr _ t -0Syr
t -0.5yr f =O. lyr \, ‘\
hydrostatic pressure
(0)
Depth
(b)
Fig. 20. Pore pressure profiles in clay layer at different times at sections f and 11(no creep, no surface load).
rock. The pore pressure changes rapidly and becomes stationaryin the crackedrock andthe sandlayer at t = 0.02 years (yr) (some days). The pore pressure profile at t = 0.02yr is shownin Fig. 16.In the clay layer very little has happened so far. In the secondphase a surface line load W = 10 kN/m’is appliedat the upper surface, Fig. 13.This loadingis made by a fictitious time procedure. Duringthis phase a local plasticzone is obtained under the load at a depth of about 2.5m.
The third phase is a consolidationprocess in the clay layer whichcauses force redist~butionin the frame. The total loadingprocedure is shown in Fig. 17. A number of cases for the leakage problem is calculated; with or without surface load, with or without creep, linear or nonlinear material. An additional calculation of where swellingis also made, the load is a pore pressure increase at the lower clay surface from i~l~ation of the same rn~~de as the decrease shown in Fii. 16.Only some of the results are discussedbelow.
592
KENNETH RUNESSONet al. 108
IO’
elastic,
no creep,
hydraulic load
elastic- plastic. no creep, hydraulic load elastic-plastic.
creep,
elastic-plastic. surface load
no creep,
elastic-plastic, increase
no creep,
settlement (m)
IOyr
hydraulic
_._..~
IO9
+~~~~~
.._
_
load
hydraulic hydraulic
load load
Fig. 21. Settlement at symmetry axis for different cases.
Plots of the displacements at two different times are shown in Fig. 18a and b. It may be observed that a yield hinge appeared in the frame at the time t = IOyr. The spread of the yield zone for the same times t = 0.5 yr and t = 10 yr is shown in Fig. 19. We observe that at t = 10yr almost the whole soil is plastic, which is a direct result of the fact that OCR is near 1 (OCR = 1.1) which means that the soil is slightly overconsolidated. Figure 20 shows the total pressure profiles in the case with no creep at different times at two sections I and II, indicated in Fig. 13. Figure 21 shows the settlement at the symmetry axis for different cases. In Fig. 21 it is observed that the case with no creep gives an asymptotic total displacement, which is not the case if creep is accounted for
7. CONCLUSION
The integrated analysis of sand-clay layer systems and the superstructure is efficiently performed according to the given theory. To obtain reliable results from the time-dependent settlement analysis the constitutive model for clay must include plasticity and creep. The adopted time-integration method gives unconditionally stable solutions and the time steps are chosen automatically based on the estimated truncation error. The numerical examples demonstrate that the used algorithm allows for the solution of complicated interaction problems, such as analysis of landslides, analysis of consolidation problems due to tunnel leakage, and analysis of swelling due to ground water injection.
REFERENCES I. B. Goshy, Soil-foundation-structureinteraction. ASCE, J. Struck Div. 104,(STS), 74s761 (1978). 2. H. L. Davidson and W. F. Chen, Nonlinear response of
drained clay to footings. Compul. Structures 8, 281-290 (1978). 3. 0. C. Zienkiewicz and P. Bettis, Fluid-structure dynamic interaction and wave forces. An introduction to numerical treatment.Int. J. Num. Mefh. Engng 14, l-16 (1968). 4. R. S. Sandhu, Fluid flow in saturated porous elastic media. Univ. of California, Berkeley, 1%8. 5. K. Runesson, On non-linear consolidation of soft clay. Chalmers Univ. of Techn., Dept. of Struct. Mech., Pub178:1, Gateborg, 1978. R. Larsson, Basic behaviour of Scandinavian soft clays. Swedish Geotch. Inst., Rep. No. 4, Linkiiping, 1977. P. Perzyna, Fundamental problems in viscoplasticity. Adu. Appl. Mech. 243-377 (1%6).
K. Terzaghi, Theoretical Soil Mechanics. Wiley, New York (1943). G. Sillfors. Preconsolidation pressure of soft, high-plastic clays. Chalmers Univ. of Techn., Geotechn. Dept., GBteborg, 1975. IO. K. Runesson and K. Axelsson, An anistropic yield criterion for clay. Proc. Inr. Conf. on Finite Elements in Nonlinear Solid and Struct. Mech. Geilo, 29 Aug.-I Sept. 1977,Tapir, Trondheim. II. N. Janbu, Foundations in Soil Mechanics (in Norwegian), Tapir, Trondheim, 1970. 12. A. Singh and J. Mitchell, General stress-strain-time functions for soils. J. Soil Mech. Found. Div., ASCE 94,214 (1%8). 13. P. Maripuu and J. Berntsson, The landslide in Tuve-a geohydrological review. J. Swedish Sot. Ciuil Engrs No. 8-9, 21-24 (1978).
14. K. Runesson, H. Tbnfors and N.-E Wiberg, GEOFEM, A Computer Program for Finite Element Analysis of Geotechnical Problems, Chalmers Univ. of Techn., Dept. of Struct. Mech., Giiteborg (in preparation). IS. A. Wolfbrandt, A study of Rosenbrock Processes with respect to order conditions and stiff stability. Chalmers Univ. of Techn., Dept. of Comp. Sci., Giiteborg, 1977. 16. K. Runesson, N-E. Wiberg and A. Wolfbrandt, A new Rosenbrock-type method applied to a finite elementdiscretized nonlinear diffusion-convection equation. The Mathematics
of
Finite
Elements
and
Applications
MAFELAP 1978.Academic Press, London (1979).
Ill,
(consolidation)or undrained conditionis considered.Transient behaviouris investigated. Yield hinge theory is adopted for the structural members while visco-plasticity, plasticity or a combination of visc&plasticity and plasticity are adopted for the soil skeleton to account for inviscid yield and creep. Adopted yield criterions are of cone type and according to the Critical State Theory. Finite element calculations for the stability analysis of a natural slope in clay and the interaction of a frame resting on a consolidating elastic-plastic subsoil are presented. 1. INTRODUCTION
Many problems include subregions with quite different physical behaviour. However, they interact and give mutual dependence, and it is of utmost importance to account for this in a unified analysis. Soil-structure interaction on structural level is described for instance by Goschy[l], and for a single structural member by Davidsson and Chen[2]. In Zienkiewicz et a/.[31 fluidstructure dynamic interaction is described. In interaction problems, substructuring technique is very convenient both for the problem description and for the numerical calculation. It is convenient that each substructure contains a special problem based on a physical law. The special problems (substructures) are coupled by interfaces by which special physical phenomena can be handled, such as crack, slide, pore pressure discontinuity etc. In this paper the interaction problems in soil due to transient flow are investigated. The paper starts with a review of some general theory for saturated soil which has to be treated as a two-phase material. The motion of the two phases (the soil skeleton and the pore water) is coupled due to interactive forces. The basic mechanical relations for the coupled problem are given elsewhere, e.g. Sandhu[4], Runesson[S]. In judging the magnitude and progress of the settlements of a soil layer, the choice of a reliable constitutive model for the soil skeleton is of outstanding importance. For soft clay it seems appropriate to include plasticity in the constitutive equations. Creep effects are important in many Swedish clays, especially when the organic contents are significant. There seems to be experimental evidence for the use of visco-plasticity as observed from undrained triaxial tests, Larsson[6]. However, though
the creep rate has been tested under fairly restrictive stress conditions, use will be made in this paper of an associated creep law according to the theory by PerzynaV]. One type of interaction is the interaction between layers of silt (or sand) and the embedding clay. Because of the larger permeability in the silt (sand) layer, the pore pressure may increase, for example due to heavy raining, much more rapidly in the silt than in the clay. As a result the effective stresses decrease first in the silt and successively in the clay. thus reducing the strength. A number of landslides in natural slopes in Sweden during recent years have made this problem paramount, Fig. la. Interaction of different types of soils claims for the use of different types of constitutive properties (e.g. yield conditions) in the same problem. In most construction problems in practice it is not sufficient to analyse the superstructure and the soil separately. In general it is more important to account for true interaction in the case of a slab-wall structure than a beam-column structure, which is more flexible. The redistribution of internal forces caused by the interaction phenomena with the soil is time dependent (consolidation and creep). Long term effects are extremely important in gravity structures. In Sweden many aquifers consist of fairly thin sand layers, which are confined from above by a thick deposit of soft clay and from below by solid rock. An important consequence of a disturbance of the hydraulic equilibrium is the occurrence of pronounced settlements at the ground surface. It is obvious that a complete analysis must include both the transient see page in the aquifer itself as well as the more delayed settlements and pore
’ settlements
clay
rock
Fig. 1. Interaction problems. 581
^ ” \ tunnel
582
KENNETHRUNESSON et al.
water flow in the clay and, finally, the behaviour of the structure under the time dependent settlements. Figure lb shows a typical situation where a frame interacting with the soil is subjected to hydrodynamic loading from leakage into a rock tunnel in the bedrock below. 2. MECHANICS FOR SATURATED SOIL
2.1 Effective stress hypothesis Saturated soil may be idealized by a mechanical model consisting of basically two phases, a solid phase (the soil skeleton) and a fluid phase (the pore water), mixed together. According to the effective stress principle by Terzaghi[8] the total stress &ii is equal to the difference between the effective stress ait in the soil skeleton and the pore water pressure p (tension stresses positive). Crij = Uij - p&j*
(1)
It is to be noticed that stresses are measured per bulk area. By definition all deformations in the soil are coupled to the effective stresses.
2.2 field equations Adopting Cartesian c~rdinates we may summarize the basic relations for saturated soil. Before any further approximations are introduced, the equations are valid both for sand and clay. ~uiiib~um equations: - (Fij,j - Uij)
t p,i = 0
(2)
meability tensor, and h,, are functions expressing the viscoplastic strain rate. Explicit time hardening has to be contained in hki for an adequate description of the behaviour of many Swedish clays. 2.3 Boundary conditions The boundary of a finite body of saturated soil is divided into two parts in two different ways. On disjunct parts relevant boundary conditions have to be imposed on displacements ui and tractions Ti as well as on the excess pore pressure p and the draining flow Q.
3. NON-LINEAR CONSTITUTION OF THE SOIL SKELETON
3.1 General Because of the great variety in ~haviour of soil depending on sedimentation environments, porosity, water contents, organic contents, cementing interactive forces, etc., there is little possibility to unify the behaviour of soil in one single constitutive model. In particular, many Swedish clays are soft with high porosity and contractant under deviatoric loading beyond the preconsolidation stresses. Even under rapid, undrained loading a dramatic yield is obtained at the state of preconsolidation, S~liforsr9]. In some cases the behaviour is almost perfectly plastic. In addition, creep effects are significant in the overconsolidated region, even for moderate additional stresses. These properties indicate that the choice of an elastic-viscop~asti~-plastic model is well motivated. Pure creep and inviscid plasticity are obtained as special cases of the more general model.
where
The weight of the skeleton represented by the initial volume forces II? is equilibrated by the initial stresses &. Additions loads are applied only on the boundary, thus making (2) homogeneous. Further, p is the excess pore pressure (diffusive pressure), since the hydrostatic pressure equilibrates the pore water bulk weight. Continuity condition: & + 4k.k = 0 (3) where rkk is the volumetric strain rate and qi are the pore water velocities relative to the soil skeleton. A dot denotes differentiation with respect to time. Equation (3) arises from the assumption of the pore water being completely incompressible. Kinematic equations (small strains): Cij= !$U,.j + Uj.i), qi = - P,i
(4)
where displacements are denoted by tci,strains by +, and diffusive forces by gi. Constitutive equations for the soil skeleton tincluding elastic-viscoplastic-plastic behaviour) and for the pore water flow (Darcy’s law): Cij = S&(&f - hki), (3 I s$, = S$&%z,,, &), hkj = h&&n, tf qi = Kigi
(6)
where S& are components of the efasti~-plastic tangential stiffness tensor, Kii are components of the per-
3.2 Elastic-viscoplastic-piastic modes The rheologicat model for the material model used is given in Fig. 2a. Associated with the two sliders are two families of closed surfaces in stress space, Fig. 2b, the inviscid yield surface (a = I) and the quasistatic yield surface (a = 2), defined by f~=O,f~=f~(c+ii,K,),(y=1,2
(7)
where K* are h~deni~ (softeni~) parameters. In the present application the inviscid yield surface f, = 0 is assumed to be hardening (softening). The plastic part of the volumetric strain is chosen as the hardening measure K$,i.e. restriction is made to pure density hardening. The quasistatic yield surface f2 = 0 is assumed to be perfectly plastic so that ~~ is irrelevant. The functions fa are similarly defined such that the surfaces f_ = 0, i = 1,2, completely coincide for a certain amount of softening in f, = 0. While f, I 0 must always hold, there is no such restriction on &. Thus, f:! = c, c > 0, define subsequent loading surfaces. Depending on the position of the stress point in stress space different characteristic stress states are obtained, Fig. 2b. The kinematic decomposition of the strain rate is (small strains) g=~;+&+&
(8)
Using associated flow rules (which may be questionable for certain choices of fa) both for the plastic part, i&, and for the creep part, i;, in accordance with Perzyna[7], we obtain the consti~tive law in Stiffness form according to (5) with the explicit expressions for
583
Interaction problems in soil structure mechanicswith materialnonlinearity 412
T
, elastic-wscoplastlc-plastic elastic-wscoplastlc
k, 2 k2 (4
(b) Fig. 2. (a) ID-rheological model, (b) loading surfaces in stress space.
the tangential stiffness (9) where
D=af,sf, afl_
-_afl afl &Yij "kra~krawk) aukk
and for the creep rate
hij
=
/I I
v(t)iO(F)>$g 1,
3.4 Hardening Based on experimental observations the material behaviour of silt and sand can be described as shear strain hardening of the extended von Mises yield criterion according to (11) with no cohesion, i.e. pb = 0. Initial yield corresponds to a cone with the slope M”’ = M., Fig. 4a. The material hardens in a drained biaxial test as in Fig. 4b, and the ultimate stage is defined as a cone with the slope M,. Intermediate cones are defined by the slope M as a function of the deviatoric plastic strain l Qp. (13)
(10)
8,
where y and @ are creep functions, F is a non-dimensional variable obtained by scaling f2 with a reference parameter. The length of a second order tensor ‘I7 is denoted ]&I = (Tk,Tk,)“*. Further, for a SC&r X, (X)=X forxrO,(x)=Oforx
(11)
containing extended von Mises (compression cone, extension cone), Drucker-Prager, etc., Fig. 3a. Surfaces associated with the Critical State (CS) theory are defined by, Fig. 3b, f~=4*tM*@,-p6)@)-p~),a=1,2
where fi(e/)=
-&)i)
(14)
with the notations ti(eqP) = 4(eq%rb M M’=&j,M,=
-
M, z- 3GeYp 1- M,/3’ 2~0 ’
The relation between the deviatoric plastic strain E,” and the deviatoric stress 4 in a drained biaxial test is shown in Fig. 4b, where the non-linear part is assumed to be a hyperbola. The test starts at the isotropic consolidation stress pb. The shear modulus in the elastic region is G. A hardening rule for clay suggested by Janbu[l I] is used in combination with yield surfaces associated with the Critical State theory, (12),
(12)
where, in particular, p)) is the isotropic preconsolidation pressure, pb is a cohesion parameter, and M is the slope of the Critical State Line (CSL), defined as h = q-M(p’-p6)=0. An anisotropic extension of (12) has been developed recently, Runesson and Axelsson[ 101,Runesson[S].
Ii& --Et X&*+2(&
PC
‘(I)_
-Pro
’ (1)
(15)
ev(kmqP)
where p :j is the initial preconsolidation pressure and k,,, is a nondimensional constant. The choice seems reasonable for Swedish soft clays. 3.5 Creep Recent tests (undrained) on Swedish clays, Larsson[6],
Q
,h=O /
,&
I
kSU , anisotropic
Fig. 3. Isotropic yield criterions: (a) Extended von Mises etc., (b) Critical State.
locus
KENNETH RUNESSON et al.
584
(a)
(b) Fig. 4. Yield criterion for silt. indicate a linear creep law in the overconsolidated region
while an exponential dependence similar to the one by Singh and Mitchell[l21 seems reliable for stresses in the post-yield region. In addition, we have to take explicit time hardening into account. Details about this can be found in Runesson[S]. 4. UNIFIED NUMERICAL ANALYSIS
4. I General A unified finite element analysis is devised for the soil-structure interaction problems discussed above. Restriction is made to ZD-cases. One very important problem is the stability of natural slopes, Fig. 5. In Sweden these problems are of current interest because of recent landslides in high sensitive plastic clays, Maripuu and Berntsson[l3], have caused very severe consequences like destroyed buildings and lost lives. The problems have a 3D-character, but in many cases it seems reasonable from the physical point of view, and necessary from the computational point of view, to approximate the problem as a 2D plane strain problem. Plasticity and creep together with the magnitude of the pore water pressure play an important role for the overall behaviour. Because of the softening character of the clay, local collapse of the clay may cause a progressive collapse. The local collapse may be caused, for example, by creep or increase of the pore water pressure. The small “local superstructure” has very little influence on the total performance of the natural slope, Fig. Sb.
It is important to analyze the influence of rapid change of the pore water pressure in the silt or sand layers embedded in the clay, Fig. 5b. Increase of pressure may be caused by a period of heavy raining. In the example in Section 6.1 a preliminary calculation based on a one phase material model is shown. Another important problem concerns the long term, time-dependent behaviour of thick clay deposits due to consolidation and creep, which cause settlements and redistribution of stresses in interacting structures. As a result severe damages on buildings may occur. Figure lb shows a situation where there is a groundwater leakage into the rock tunnel from the sand layer. The resulting pore pressure drop in the sand layer will, however, very soon become stationary. This decrease of pore pressure in the sand layer starts a consolidation process in the clay layer, which results in settlements at the ground surface. The frame interacts with the soil so that the settlements cause additional displacements and stresses in the frame. A third case is to analyse the swelling problem. In many cities a lowering of the groundwater table has led to consolidation and thus settlements. By infiltration it is possible not only to stop the consolidation but also induce the surface to raise again to a certain extent. The calculation is valid only for a moderate increase in the pore pressure. The problem is nonlinear and time dependent. The universal parameter is the time t. Loads W, boundary conditions and constitutive parameters may be time dependent, Fig. 6a.
i Lic+Ls
Fig. 5. (a) 3D-natural slope (valley), (b) sdt layer in natural slope, ZD-plane strain.
W
w,
‘,
_____--
I
(a)
\ W,i/i, Fig. 6. Real- and fictitious time.
(b)
585
Interaction problemsin soil structure mechanics with material nonlinearity In order to handle an instantaneous load a fictitious time parameter i is defined. The instantaneous load is applied to the structure as W = W,di, at the real time t,, Fig. 6b. Thus, the action spectra can be applied to the structure in the time space (t, t), see Fig. 6c. 4.2 FE-discretization in space
A discretization in space of an appropriate variational formulation leads to the semi-discrete Galerkin method for problems with time dependence. A system of ordinary differential equations in time is obtained. The space discretization used in the computer program GEOFEM [ 141,for calculation of the 2D-problems shown below, is based on a parametric finite element with a biquadratic displacement approximation and a bilinear excess pore pressure approximation in local coordinates. The element geometry is mapped from quadratic regions in the local coordinates by biquadratic functions. For coupled consolidation such an element may be termed quasi-isoparametric, while for the pseudo-theory (ground water flow analysis) the element will be superparametric.
terms have the required regularity to match the theoretical order of accuracy. As the coefficient matrix will be constant during a time-interval, an extra stage in the method corresponding to the estimation of the truncation error only requires an extra back-substitution provided an elimination solution method is used, which is the case in the computer program GEOFEM. The particular method adopted is called (2,4)-W. It is a 2-stage method of order 2. Applied to the problem
where f is a column matrix of the same dimension as j, the method is defined by j,=&&3k*) W(& 0, P)g, =f(&> to) W(At, (Y,P)k,=f($,t;AtE,,
4.3 Numerical time integration The non-linear initial value problem arising from the FE-discretization of the consolidation problem can be formulated in matrix form as _d6 A@&=-B6t
IC(O)=co.
V(t)tA(& t)
WI
The nodal values of the displacement S and the excess pores pressure p’ are collected in the matrix 5, while the applied load U diffecentiated with respect to time and the known drainage P are collected in V. A pseudoload due to creep, f?‘(tj, t), is contained in A.
The square matrix A, which is state dependent, and B, which is constant, are defined by
t,,t;At
(18) -;AfaPk,
W(At, (Y,P)=Z-AtaP where P is an arbitrary real quadratic matrix and a is the inverse of the largest root to the second Laguerre-polynomial, i.e. a = 1 - l/6(2). The method has a built-in local error estimate of third order. Details may be found in[15]. Due to the certain kind of non-linearity introduced in (16) at the adoption of an elastic-viscoplastic-plastic model, some modification is made at the practical application to the problem of consolidation. This is motivated from a solution economy point of view. The main modification is that the tangent stiffness matrix S is evaluated for t& only and is kept unchanged through all stages during a single time step. Analogeously to solid analysis, the fulfilment of the yield condition in the end of the time interval results in a lack of equilibrium, which is compensated for by applying a fictitious residual load. This load gives rise to additional displacements and pressures. An iterative procedure leads to the simultaneous satisfaction of the yield condition and the variational equilibrium equation. 5COMPUTER
PROGRAM GEOFEM
General A computer program GEOFEM [ 141for the solution of non-linear and time-dependent geomechanical problems is under development. Main characteristics are: 0 Simplicity and generality in the program structure 0 New types of problems can be added in a simple way 0 Input data is given free format introduced by key words 0 The calculations are steered by simple input commands 0 Substructuring technique is used 5.1
where S is the tangential stiffness matrix for the skeleton and C is a rectangular coupling matrix. Further, K is the permeability matrix and GS is associated with a convective flow boundary condition. In the program GEOFEM a single-step method of linear semiexplicit Runge-Kutta type is used to trace the solution in time. The method used is basically a so called W-method, a class of methods which has recently been proposed by Wolfbrandt [151. W-methods are convenient for stiff systems like (16). An automatic prediction of the next time step is based on the calculation of the principal error term, which is determined by a comparison of the results obtained from two methods of different order (imbedding technique). Hence, for small time steps a reliable estimate for the local truncation error is possible, provided the non-linear
5.2 Substructuring. Interfaces The structure is described by a number of substructures coupled by interfaces. The interfaces may either be a logical coupling, or a coupling containing stiffness (or other material) properties. By use of the interfaces cracks, friction, sliding, discontinuities in pore pressure,
KENNETH RUNPSONel 01.
586
and so on can be described. Each substructure has a simple description and the output is obtained separately for each substructure. The substructures may have different mesh sizes, which means that master nodes are used at the interfaces in order to fit the different meshes together. Two cases of interaction from soil-structure mechanics are shown in Fig. 7a and b. Master variables li” in the master nodes are introduced by the transformation (ti refers to a certain substructure)
The normal stresses in situ in the xyz-system (~~~,u~~,u~~) are assumed to be the principal stresses. The total vertical stresses and pore pressures are calculated as follows, Fig. 9:
a=o= i p. =
g=fi*“*
(19)
The stiffness relation for interface nodes are obtained by changing the stiffness relation 5% = U to D’SDC” = o”, Ljm=Dto.
(20)
These multiplications are made in condensed form. Addition of stiffness properties in the interfaces yield the final equations
[
sl;,Ys2;s;][ ;;:I=[;;:]
(21)
where index 1 and 2 refers to substructure 1 and 2, and I refers to the interface. 6.~CAL 6.1.Stability analysis
EXULT
6.1.1 Problem. Basic data. A stability analysis of a natural slope from Tuve in GGteborg, Sweden[l3], is performed. The stope is shown in Fig. 8. We want to estimate the safety factor for the slope. The safety factor is here defined as the factor fi, which multiplied with the total weight of the slope gives instability. The influence of a silt layer on the safety factor is studied.
2 E LZ&zrl z E [z,, &I z E [Z%z,l
P&z -Lo), PO& - zo)+ PC?& - z,), PO&, - zo)+Plzg(&-~I),
pgtz - z2), .?E [zo, ZJ apug@ - z,) + APO,p,,g = 10kNlm’, z E IZI, zzl z E I&, ZJ. i %vg(z - ZI),
The pore pressure p. is dependent of the flow conditions, because of an inclining ground water level. This is accounted for by the factor a. For pure hydrostatic conditions a = 1. Figure 9a shows the distribution of the total vertical stress aVo and the pore pressure p. at different depths, where Ape means a drastic increase of pore pressure in the silt layer (which may be caused by season variations). The local increase of the pore pressure in the silt layer may be approximately calculated by simple transient aquifer analysis. This flow problem becomes stationary in a short time (some weeks or a couple of months). The failure will not take place in the silt layer but in the clay nearby, which justifies the calculations made here. The effective stresses in situ will thus be cr:o=a;o-PO,
(22)
(aLo= -aLo)
0:o = KoaLo, (&o= -Cr:o= -&o). 6.1.2 Used method. The analysis is here performed as an undrained calculation for a one phase material (total stress analysis). For an isotropic elastic material the effective mean pressure p’ is then constant during the / frame
f ting
sltt layer master node (Interface with physical
la)
(b)
Fig. 7. Interfaceswith master nodes.
r
+...__._
-_,x
i z
12m
----_.-._-..~-.____
250m
__
I
Fig. 8. Natural slope
in
Tuve.
4
Interaction
problems in soil-structure mechanics with material nonlinearity
587 kN/m
50
I,(X)
z,(x)
I
1
”
”
,,
@‘“=I5
1
t
9.2
i
4
P0,g=17.5kN/m2 I
z(m)
04
(4
Fig. 9. Distribution of: (a) vertical stresses Q, and pore pressure p0 and, (b) undrained cohesion c, for a vertical profile. loading. The strength can be represented with a deviatoric parameter which is a function of the effective mean stress in situ. Von Mises yield condition is chosen which is a special case of the extended von Mises (DruckerPrager) with
responds to the true deviator strength q. = 2c, for the initial stress state, is determined from (25) by setting @, = 0, Fig. IO. The condition qd = q. for p’ = PAgives, if sin & < 1, c,, = pb sin & t cd cos &,
c$,,= 0.
(26)
C#J = r#J”= 0, c = c,. The determination of the undrained strength c, is demonstrated below. For synchronic loading the proportionality factor serves as a fictitious time variable. The nonlinear initial value problem from the FE-descritization can be formulated in matrix form as (a part of (16))
The effective mean stress pb in situ is defined by, the state stress (TX0 = uyo = (aLo, o:o, uL0) with K do, ( do = - dd as p;=
_1/ it/c = 3
;(I+2K&r:o.
(27)
For clay it is also possible to calculate the C.-value from the undrained shear strength 7fU obtained from (23) vane borer tests. with the constraints of weak equilibrium
&(a - u,,) = Uf, b, =
I ”
cu =
A&‘o d V
(24)
and the yield condition. 6.1.3 Calculation of c,. The average cone yield condition (concerning effective stresses) can be formulated as 18sin&, , 18co&, (25) qd = 9_sinz,$dp ‘m ‘d. The undrained
cohesion parameter
c., which cor-
;[(1-
Ko)*a:a + 3&"*.
It is then assumed that the normal stresses in situ remain unchanged during the test. For the problem in Fig. 8 with the ground water table at the ground surface, the data according to Fig. lob have been chosen. Inserted in (26) it gives the undrained cohesion c. shown in Fig. 9b. A second calculation is made for an increase Appoof the pore pressure, which gives the change in the c.-value AC. = f( 1t 2K,,)( - Ap,,)sin c#+, 6.1.4 Numerical
calculations.
Figure lla shows the
4
T
ApI
Clay:C,,=IOkN/m’ Ko= 06 &= 20” a = I.1
Sand:
(28)
cd=0 K,J=0.6 &=36’ Y 50.495 (b)
P’
(a)
Fig. IO. (a) Simulation of yield condition with undrained parameters, (b) soil data.
KENNETH RUNESSON ef al.
/?=I.55
1.75
1.55
2.15
1.95
2 15
(b)
(a) Fig. 11. Slope: (a) displacement picture for p = 1.65(displacements magnified I5 times), (b) gowth of plastic zone.
finite element mesh and the displacement picture for the load factor p = 1.65. Instability is obtained for the load factor 2.15. The growth of the plastic zone is shown in Fig. lib. Figure 12 shows the displacement picture and growth of the plastic zone for the case with a drastic pore pressure increase of 40 kN/m* in the silt layer giving AC. = - 10kN/m*. Instability in this case is found for the load factor 1.7. In particular, we observe the large shear deformations in the silt layer and the development of the “circular” yield zone. A local zone shown in Fig. 12 (as a shadowed area) with a low shear strength value diminishes the load factor to 1.4. We observe that the failure load is very sensitive to local areas with low shear strength. A more sophisticated and true analysis ought to include creep and softening behaviour of the soil.
6.2 Soil-frame stress interaction and sand-clay fZow interaction The leakage problem in Fig. 13 is studied under consideration of different material properties and different loadings. Both the flow interaction between the sand and the clay layer and the stress interaction between the soil and the frame have to be accounted for. The long term effects from consolidation and creep are studied. The soil displacements and the displacements and redistributed forces in the frame are obtained as results. 6.2.1 Basic data. The soil profile in Fig. 14a is considered under two different loading alternatives. The initial pore pressure distribution is assumed to be hydrostatic. The vertical effective pressure (compression stress) in
Fig. 12. (a) Displacement picture for p= 1.65(displacement magnified 15 times), (b) growth of the plastic zone.
589
Interaction problems in soil-structure mechanics with material nonlinearity
sym
1
7m
*
frame,
c
Sand
o .
- k=lO-s
.
II
I
7m
I
T
I I
I I
.,* II
.
~
III.
o
. .
.
3m
Ah,hhhh
IOm
3m
3m
Fig. 13. Clay-sand flow interaction and frame-soil interaction.
IOm Soft
clay
Fig. 14. (a) Soil profile, (b) initial stress a:,, and preconsolidation stress 0: distribution with depth in situ is assumed as a:o=&,
p’ = Pm - PW.
A normal value for Swedish clay is p,,,g=16kN/m3 giving p’g =6 kN/m’. It is assumed that the preconsolidation pressure ail), which is the equivalent to ~21)for the anisotropic extension of (12),[5], can be estimated as -‘(I) = a:( 1+ 2&)/3 PC where, in the present example, Z&=0.8 is chosen. The distribution of a:, which is the vertical preconsolidation stress obtained at an oedometer test, is assumed as, Fig. 14b, I 2.0~~kN/m*, 01z12m (Tc=I za kN/m*, 2~z~lOm where cx is a constant. The constant value of (T: in the upper part of the clay layer corresponds to a dry crust. The quasi-static preconsolidation pressure & is assumed to be equal to p& the mean effective pressure in situ. It can also be calculated as B:“‘/OCR
where the in situ value of the overconsolidation ratio OCR is OCR+= I
2.OiYlp’gz, 0 5 z I2m ( fflp’g, 2 5 z I 1Om
Beneath the dry crust a slight degree of overconsolidation is assumed by the choice of cr=6.6 (OCR= 1.1). The elastic constants used are E=3000 kN/m’ and v=O.30. The internal friction angle is a=30 (giving M) and a small cohesion is assumed, c=S kN/m* (giving pb). To define the hardening along the normal consolidation line at &-consolidation, the value k,,,=4 has been chosen. The isotropic permeabiiities for the clay, sand, rock and cracked rock are found in Fig. 13. 6.2.2 Numerical calculations. The finite element mesh is shown in Fig. 15. The figure also shows the different types of elements; beam elements (3 variables at each node), consolidation elements (displacements li and pore pressure @), and ground water flow elements (pore pressure 6). The drilling of a rock tunnel initiates a ground water leakage from the sand layer through the zone of cracked
KENNETH RUNESSON et al.
SnuoP%sY”re
I
, soll-structure stress interaction
4 ’
r-_A clay layer
I id -~- ~1
slow transient flow
11
, sand-clay flow interactlon
! P--i w 4 sand and rock, rapid + c transient flow
I I
j
I
1 +I3
___ It-L Fig. 15. Finite element mesh. s
n
clay
t--o
. .. sand ,. . hAx
pore pressure proflle stationary solutlon for the decrease ,of the pore pressure at the upper vurface of the sand layer pore pressure 1 ,
rock symmetry 0x15 after 0.02years
\
hydrostatic
\ \J
/
-
\
: _t,t_+_. 100
3
\
:
:
I
pressure
:
:
,pkNh2
4
300
200
100
tunnel
p kN/m’
Fig. 16. Pore pressure protiles at t = 0.02yr (no creep, no surface load). W
,t tunnel drilling
’ mainly flow in sand layer
consolidation external loadmg
Fig. 17. Loading procedure in real and tictitious time.
591
Interaction problems in soil-structure mechanics with material nonlinearity
displacements O.lJm
Fig. 18. Displacements at (a) t = 0.5 yr and (b) t = 10yr (no creep, no surface load). t=O.Oiyr
sym
I
Fii. 19. Spread of yield zone at t = 0.5 yr and t = 10yr (no creep, no surface load). Section 20
40
Section II -11
I-I
60
60
too :
, p(kN/mz)
\.
k
,=aa
6
Depth (ml
t = IOyr
=IOyr _ t -0Syr
t -0.5yr f =O. lyr \, ‘\
hydrostatic pressure
(0)
Depth
(b)
Fig. 20. Pore pressure profiles in clay layer at different times at sections f and 11(no creep, no surface load).
rock. The pore pressure changes rapidly and becomes stationaryin the crackedrock andthe sandlayer at t = 0.02 years (yr) (some days). The pore pressure profile at t = 0.02yr is shownin Fig. 16.In the clay layer very little has happened so far. In the secondphase a surface line load W = 10 kN/m’is appliedat the upper surface, Fig. 13.This loadingis made by a fictitious time procedure. Duringthis phase a local plasticzone is obtained under the load at a depth of about 2.5m.
The third phase is a consolidationprocess in the clay layer whichcauses force redist~butionin the frame. The total loadingprocedure is shown in Fig. 17. A number of cases for the leakage problem is calculated; with or without surface load, with or without creep, linear or nonlinear material. An additional calculation of where swellingis also made, the load is a pore pressure increase at the lower clay surface from i~l~ation of the same rn~~de as the decrease shown in Fii. 16.Only some of the results are discussedbelow.
592
KENNETH RUNESSONet al. 108
IO’
elastic,
no creep,
hydraulic load
elastic- plastic. no creep, hydraulic load elastic-plastic.
creep,
elastic-plastic. surface load
no creep,
elastic-plastic, increase
no creep,
settlement (m)
IOyr
hydraulic
_._..~
IO9
+~~~~~
.._
_
load
hydraulic hydraulic
load load
Fig. 21. Settlement at symmetry axis for different cases.
Plots of the displacements at two different times are shown in Fig. 18a and b. It may be observed that a yield hinge appeared in the frame at the time t = IOyr. The spread of the yield zone for the same times t = 0.5 yr and t = 10 yr is shown in Fig. 19. We observe that at t = 10yr almost the whole soil is plastic, which is a direct result of the fact that OCR is near 1 (OCR = 1.1) which means that the soil is slightly overconsolidated. Figure 20 shows the total pressure profiles in the case with no creep at different times at two sections I and II, indicated in Fig. 13. Figure 21 shows the settlement at the symmetry axis for different cases. In Fig. 21 it is observed that the case with no creep gives an asymptotic total displacement, which is not the case if creep is accounted for
7. CONCLUSION
The integrated analysis of sand-clay layer systems and the superstructure is efficiently performed according to the given theory. To obtain reliable results from the time-dependent settlement analysis the constitutive model for clay must include plasticity and creep. The adopted time-integration method gives unconditionally stable solutions and the time steps are chosen automatically based on the estimated truncation error. The numerical examples demonstrate that the used algorithm allows for the solution of complicated interaction problems, such as analysis of landslides, analysis of consolidation problems due to tunnel leakage, and analysis of swelling due to ground water injection.
REFERENCES I. B. Goshy, Soil-foundation-structureinteraction. ASCE, J. Struck Div. 104,(STS), 74s761 (1978). 2. H. L. Davidson and W. F. Chen, Nonlinear response of
drained clay to footings. Compul. Structures 8, 281-290 (1978). 3. 0. C. Zienkiewicz and P. Bettis, Fluid-structure dynamic interaction and wave forces. An introduction to numerical treatment.Int. J. Num. Mefh. Engng 14, l-16 (1968). 4. R. S. Sandhu, Fluid flow in saturated porous elastic media. Univ. of California, Berkeley, 1%8. 5. K. Runesson, On non-linear consolidation of soft clay. Chalmers Univ. of Techn., Dept. of Struct. Mech., Pub178:1, Gateborg, 1978. R. Larsson, Basic behaviour of Scandinavian soft clays. Swedish Geotch. Inst., Rep. No. 4, Linkiiping, 1977. P. Perzyna, Fundamental problems in viscoplasticity. Adu. Appl. Mech. 243-377 (1%6).
K. Terzaghi, Theoretical Soil Mechanics. Wiley, New York (1943). G. Sillfors. Preconsolidation pressure of soft, high-plastic clays. Chalmers Univ. of Techn., Geotechn. Dept., GBteborg, 1975. IO. K. Runesson and K. Axelsson, An anistropic yield criterion for clay. Proc. Inr. Conf. on Finite Elements in Nonlinear Solid and Struct. Mech. Geilo, 29 Aug.-I Sept. 1977,Tapir, Trondheim. II. N. Janbu, Foundations in Soil Mechanics (in Norwegian), Tapir, Trondheim, 1970. 12. A. Singh and J. Mitchell, General stress-strain-time functions for soils. J. Soil Mech. Found. Div., ASCE 94,214 (1%8). 13. P. Maripuu and J. Berntsson, The landslide in Tuve-a geohydrological review. J. Swedish Sot. Ciuil Engrs No. 8-9, 21-24 (1978).
14. K. Runesson, H. Tbnfors and N.-E Wiberg, GEOFEM, A Computer Program for Finite Element Analysis of Geotechnical Problems, Chalmers Univ. of Techn., Dept. of Struct. Mech., Giiteborg (in preparation). IS. A. Wolfbrandt, A study of Rosenbrock Processes with respect to order conditions and stiff stability. Chalmers Univ. of Techn., Dept. of Comp. Sci., Giiteborg, 1977. 16. K. Runesson, N-E. Wiberg and A. Wolfbrandt, A new Rosenbrock-type method applied to a finite elementdiscretized nonlinear diffusion-convection equation. The Mathematics
of
Finite
Elements
and
Applications
MAFELAP 1978.Academic Press, London (1979).
Ill,