- Email: [email protected]
sandy soil with a concrete lining
Engineering Structures 22 (2000) 707–722 www.elsevier.com/locate/engstruct
Nonlinear analysis of tunnels in clayey/sandy soil with a concrete lining S.A.S. Youakim a, S.E.E. El-Metewally b, W.F. Chen b
c,*
a Arab Consulting Engineers, Cairo, Egypt Structural Engineering Department, El-Mansoura University, El-Mansoura, Egypt c School of Civil Engineering, Purdue University, West Lafayette, IN 47907, USA
Received 22 January 1998; received in revised form 12 January 1999; accepted 12 January 1999
Abstract In the nonlinear analysis of tunnels embedded in clayey or sandy soil with concrete lining, nonlinearity arises due to three parameters. The first parameter is the nonlinear material behavior of the lining and the surrounding soil. The second parameter is a result of the large deformations effect. The third parameter is the nonlinear nature of the contact between the concrete lining and the surrounding soil (essentially the slippage between the two materials). In this paper the three parameters of nonlinearity are modeled and implemented in a finite element program. 2000 Elsevier Science Ltd. All rights reserved. Keywords: Contact problems; Large displacements; Soil plasticity; Soil–structure interaction; Tunnels; Cracks
1. Introduction The response of an engineering problem such as soil– structure interaction subjected to static and/or dynamic loading is influenced significantly by the existence and behavior of the contact between the two different materials. Over the past years, analysis of soil–structure interaction problems has been carried out by assuming a complete bond at the interface, i.e., no relative motions are allowed. Although this assumption simplifies the analysis and design, it can introduce errors in the prediction of stresses and deformations. Moreover, in most of the latter researches in this field, a small deformation analysis has been assumed. But it was found that this assumption is not often true especially in the case of weak materials such as clayey soils. Another important aspect in the study of interaction problems is the nonlinear behavior of the foundation material and the surrounding soil mass. Although the application of classical limit analysis techniques to geotechnical problems gives rather satisfactory results, in many problems involving complex geometry or conditions in which knowledge of
* Corresponding author. Tel.: ⫹ 1-765-494-2254; fax: ⫹ 1-765496-1105.
deformations are important, the implementation of the advanced plasticity models may be very important. As a result of the rapid development of digital electronic computers with powerful numerical methods such as the finite element method, it is now possible to implement all the previously mentioned factors of nonlinearity in computer programs and perform a nonlinear finite element analysis for any selected problems and obtain results with a convenient degree of accuracy. The objective of this paper is to implement the nonlinear constitutive modeling of soil and concrete in addition to the assumption of large deformation analysis to a selected circular tunnel embedded in soft clayey soil or sandy soil and study how these factors may alter the predicted surface settlement or may change the stress distribution around the tunnel opening. Another intention of this paper is to study the phenomena of soil–structure interaction and how this may release or intensify the stresses in the tunnel lining. In this study, the finite element method is applied to a circular tunnel problem excavated by the shield method far away from the tunnel inlet and outlet faces. So, it is sufficient to consider the case of plane strain neglecting the advance velocity of the shield machine as it moves along the tunnel axis.
0141-0296/00/$ - see front matter 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 1 - 0 2 9 6 ( 9 9 ) 0 0 0 0 8 - 5
708
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
2. Nonlinear constitutive relations 2.1. Material modeling There exists a large variety of models which have been proposed in recent years to characterize the stress– strain and failure behavior of concrete and soil. All these models have certain inherent advantages and limitations which depend to a large degree on their particular application. In this paper the concrete material is modeled by a hypoelastic model whereas the soil is modeled by the cap model. It must be noted that in this paper our discussion is only limited to the time and temperature independent stress–strain constitutive relations. 2.1.1. Hypoelastic failure model for concrete Since concrete is an irreversible and load-path dependent material, it is appropriate to use the incremental stress–strain formulation, known as the hypoelastic model, in which the loading history effect is incorporated; thus, the load path-dependency behavior including the unloading behavior is represented. Bathe and Ramaswamy [5] developed a hypoelastic model for concrete based on a uniaxial stress–strain relation that is generalized to take biaxial and triaxial stress conditions into account. Tensile cracking and compression crushing conditions are identified using failure surfaces. The use of tensile and compression fail-
Fig. 1.
ure criteria including strain softening conditions prevents the unrealistically large stress and strain conditions predicted as in the case when using some plasticity models. With regard to compression failure, the failure envelope shown in Fig. 1 can be used. The shape of this failure surface is largely based on the experimental results by Kupfer et al. [11] after being extended to the three dimensional case. The compression failure envelope requires 24 discrete stress values to be input. Firstly, the values ip1/˜ c are input. These values define at what stress magnitudes tp1, the discrete two-dimensional failure envelopes for additional stress tp2 and tp3 are input. Then, the failure envelopes are defined by the failure stress value ijp3/˜ c (i ⫽ 1, . . ., 6; j ⫽ 1, 2, 3) that correspond to the stress magnitude tp2 ⫽ tp1, tp2 ⫽ tp3 ( is a constant), and tp2 ⫽ tp3, where tp1, t p2 and tp3 are the principal stresses at time t, with t p3 ⱕ tp2 ⱕ tp1. Tensile failure occurs if the tensile stress in a principal stress direction exceeds the tensile failure stress as expressed in Fig. 2. In this case, it is assumed that a plane of failure develops perpendicular to the principal stress direction. The effect of this material failure is that the normal and shear stiffnesses across the plane of failure are reduced, and the corresponding normal stress is released. Interested readers may refer to [5] for detailed derivation and development of the proposed model.
Triaxial compression failure envelope model.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
709
Fig. 2. Triaxial tensile failure envelope model.
2.1.2. Cap model for clayey soils The isotropic work-hardening plasticity cap model was first suggested by Drucker et al. cited in [9]. They introduced a spherical end cap to the Drucker–Prager model, in order to control the dilatancy of soil. DiMaggio and Sandler [8] proposed the concept of a strain hardening elliptical cap fitted to the perfectly-plastic (failure) surface of Drucker–Prager. An associated flow rule was employed and the movement of the cap is controlled by the increase or decrease of the plastic volumetric strain. From a theoretical point of view, the cap model is particularly appropriate to soil behavior, because it is capable of treating the conditions of stress history, stress
Fig. 3.
path dependency, dilatancy, and the effect of the intermediate principal stress. The schematic shape of this model is shown in Fig. 3. The loading function consists of two parts: 1. An ultimate failure envelope of the Drucker–Prager type which serves to limit the maximum shear stresses in the material and has the form: f 1 ⫽ ␣¯I1 ⫹ J1/2 2 ⫺ k ⫽ 0
(1)
where ¯I1 is the first invariant of the effective stress tensor; J2 is the second invariant of the deviatoric stress tensor; and ␣, k, are material constants related to the friction and cohesion of the soil.
Cap model; (a) I¯1 ⫺ J1/2 plane; (b) deviatoric plane. 2
710
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
2. A strain-hardening elliptic cap which has the form: f 2 ⫽ (I¯1 ⫺ L)2 ⫹ R2J2 ⫺ (X ⫺ L)2 ⫽ 0
(2)
where R is the ratio of the major to minor axis of the elliptic cap; and X and L define the ¯I1 value at the intersection of the elliptic cap with the ¯I1 axis and the failure function, respectively. The cap position, X, is controlled using the hardening function by DiMaggio and Sandler [8] as follows: X⫽⫺
冉 冊
1 ⑀pv ln 1 ⫺ D W
(3)
where W defines the maximum plastic volumetric compaction that the material can experience under hydrostatic loading (Fig. 4); D is a curve fitting parameter; and ⑀pv is the plastic volumetric strain. Detailed description of the cap model and the input parameters may be found elsewhere as in [6,9,10,12]. 2.2. Large deformation analysis In the present paper, the analyses are performed using Total-Lagrangian Formulation assuming both material and geometric nonlinearity considering large displacements, large rotations, but small strains. Bathe [3] gives detailed derivations of the governing finite element equations, stress–strain relations and the corresponding strain–displacement transformation matrices for different types of nonlinear analyses.
Fig. 4.
3. Constraint equations approach for soil–structure interaction
The response of a soil–structure interaction problem involving contacts or interfaces and subjected to static and/or dynamic loading is influenced significantly by the existence and behavior of the contacts. Many analyses of structure–foundation systems are carried out by assuming complete bonding at the interface; i.e., there is no relative motion. Although this assumption simplifies the analysis and design, it can introduce errors in the prediction of stresses and deformations. The main reasons for such errors are based on the fact that soil is neither homogeneous nor elastic and invariably much weaker than the structure it has to support. Moreover, experiments have shown that sliding between concrete/soil and steel/soil occurs at levels of interface shear which are significantly less than the shear limit of soil. Thus, analyses which assume a perfectly bonded interface condition over predict the shear transfer and depending on specific application, either over or under estimate structural response. Various methods have been proposed in the past in order to model the discontinuous behavior at the interface. Some workers have suggested the connection of the elements of soil and concrete to each other by discrete springs. Other applications have treated the interface by using elements of zero thickness or by thin-layered elements of small thickness. Some other methods use the approach of constraint equations or the total potential. This latter approach is adopted in this study and its main features are outlined in the following subsections.
Relationship between mean normal stress and volumetric strain.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
3.1. Concept
711
this case (tt is applied to equal dtn).
This solution approach is applicable to the analysis of contact problems involving sticking, frictional sliding and separation under large deformations. It was first presented by Bathe and Chaudhary [4]. It contains the following major aspects: 1. The total potential of the contact forces is included in the variational formulation to enforce the geometric compatibilities along the contact surfaces. 2. In the region of contact, surface tractions are evaluated from the externally applied forces. 3. The surface tractions between nodal points are employed to decide whether a nodal point is sticking, a sliding contact or releasing. 4. The number of equations due to the contact conditions is dynamically adapted to solve two equations for each node in contact if the node is in the sticking condition, and one equation if the node is in the sliding contact condition. This approach has the ability to solve contact problems whilst at the same time including large deformation effects, repeated contact and separation between the bodies at the interface. Hence, it will be considered in the present study. The approach considers two generic bodies denoted as contactor and target. In the finite element solution, the contactor contains the finite element boundary nodes that come into contact with the target segments or nodes. The basic conditions of contact along the contact surfaces are that no materials overlap can occur, and as a result, contact forces are developed and act along the regions of contact between the target and contactor. These forces are equal and opposite. The normal tractions can only excert compressive action, and the tangential tractions satisfy Mohr Coulomb’s law of friction. 3.2. Contact criteria Let tt ⫽ the developed tangential tractions along contact surfaces; tn ⫽ the compressive normal tractions; s ⫽ static coefficient of friction; and d ⫽ dynamic coefficient of friction (d ⱕ s). sticking contact:
The two bodies are assumed to be in sticking mode if 兩tt兩 ⱕ stn, which means that the frictional resistance during contact is sufficient to prevent sliding.
sliding contact:
sliding occurs between target and contactor bodies when 兩tt兩 > stn; i.e., the frictional capacity during contact is smaller than the applied tangential traction and
Tension contact:
The two bodies are assumed to separate from each other when tn ⬍ 0, and in this case both normal and/or tangential tractions are set to be zero.
4. Computer program—ADINA-96 ADINA-96 (Automatic Dynamic Incremental Nonlinear Analysis) is a computer program for the static and dynamic displacement and stress analysis of solids, structures and fluid–structure systems. The program has been written in FORTRAN 77, and can be employed to perform linear or nonlinear analysis in either plane strain, plane stress, axi-symmetric and general threedimensional analysis. ADINA-96 is written in a modular form so that any additional material models can be implemented without major programming effort. The program ADINA is the main module of the ADINA system which consists currently of the program ADINA for displacement and stress analysis, ADINAT for analysis of heat transfer and field problems, ADINA-IN for preparation and display of the input data and ADINA-PLOT for the display of the calculated solution results. The ADINA program was first developed by Professor K.J. Bathe [1]. Since then, it was subjected to various improvements in order to widen its element library, material models, and solution capabilities. The current version of the ADINA-96 program has two additional improvements added to the ADINA-84 computer program [2]. The first improvement is the implementation of the elliptic cap model for soils [10,12]. The elliptic cap model for soils was implemented after some modifications to justify the general control variables in the element library of ADINA-84 and the general analysis procedures. The second additional feature in ADINA-96 is the preprocessor program AUTOGEN [2] to assist in the preparation of the input file to ADINA by generating the model coordinates, two-dimensional element data, boundary conditions, and applied loads which balance the initial stresses in soil. In nonlinear analysis, ADINA-96 performs an incremental load-displacement analysis. After each increment of load, the computed displacement field is modified using an iterative procedure until an equilibrium configuration is reached.
5. Application to plane strain tunnel problem Soil response to construction procedures of tunnels have always been an area of great interest to engineers
712
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
since care must be taken to insure that structures adjacent to, or above, tunnel construction are not damaged by excessive ground movements. During tunnel excavation, the state of stresses in the ground continuum are significantly altered specially in soft soils. So, an immediate support must be provided to the surrounding ground; otherwise, the tunnel will collapse. On the other hand, prediction of the behavior of the tunnel lining during construction is a complicated problem due to the fact that neither the lining nor the ground around a tunnel behaves in an elastic manner. Another important factor is the nonlinear nature of the contact between the concrete lining and the surrounding soil mass. In this section, a parametric study is made to a selected circular tunnel problem concerning the type of nonlinear analysis, tunnel depth, tunnel thickness, interface behavior, lining modeling, and surcharge loads. 5.1. Input parameters for material models 5.1.1. Soil cap-model parameters A normally consolidated soft clayey soil, called the Boston-Blue Clay, is chosen to represent the soil around the tunnel cavity and its properties are presented in [6]. Humphrey and Holtz [10] used this type of soil in a nonlinear finite element study of reinforced embankments. The values of the different parameters of the capmodel for this type of soil can be found in [10,12]. 5.1.2. Concrete model parameters The behavior of the concrete lining of the tunnel is modeled by the hypoelastic model of Bathe and Ramaswamy [5] as discussed in Section 2.1.1. Table 1 summarizes the material and control parameters of the concrete model while Table 2 illustrates the parameters used to define the triaxial compression failure envelope, where E˜0, , ˜ t, ˜ c, e˜c, ˜ u and e˜u are the initial Young’s modulus, Poisson’s ratio, tensile stress, ultimate compressive stress, crushing stress, strain corresponding to crushing stress, ultimate compressive stress and ultimate compressive strain, respectively, and the tilde (~) stands for uniaxial condition.
Table 1 Material and control parameters for the concrete model Parameter
Value
E˜0 ˜ t ˜ c e˜c ˜ u e˜u
28 ⫻ 105 t/m2 (28.0 GPa) 0.18 400.0 t/m2 (4.0 MPa) 4000.0 t/m2 (40.0 MPa) 0.002 3500.0 t/m2 (35.0 MPa) 0.003
Table 2 Triaxial compression failure curves Principal stress ratios Curve number
p1/˜ c p3 at p2 ⫽ p1 p3 at p2 ⫽ 0.75 p3 p3 at p2 ⫽ p3
1
2
3
4
5
6
0.0 1.0 1.3 1.25
0.25 1.4 1.5 1.45
0.5 1.7 2.0 1.95
0.75 2.2 2.3 2.25
1.0 2.5 2.7 2.65
1.2 2.8 3.2 3.15
5.2. Problem geometry and sequence of analysis The chosen circular tunnel for this study is of 9.0 m outer diameter. This diameter is constant throughout the parametric study. In the case of a tunnel cavity without lining, a smaller diameter of 5.0 m is chosen in order to guarantee the stability of the opening, otherwise very large deformations would occur which may lead to a complete failure of the problem. Symmetry is assumed about the vertical axis so that one half the problem is considered in the analysis. The extent of the finite element mesh from the axis of symmetry is taken to be six times the diameter. This was found adequate to represent the ground subsidence because the deformations beyond this distance become nearly negligible. In the present analysis, it is assumed that the tunnel is excavated by the shield method and is located far from both the tunnel inlets and outlets, so that the assumption of plane-strain problem will be adequate. The soil continuum is represented by two dimensional 4–5 node quadrilateral elements. The concrete lining is represented by two dimensional 6-node quadrilateral elements. The contact between the lining and the soil is represented along the outer diameter of the tunnel. The nodes belonging to the soil and lining are initially assigned the same coordinates. Fig. 5 shows the whole geometry of the problem, while Fig. 6 demonstrates the finite element mesh at contact surfaces. The sequence of the analysis is organized as follows: first the problem is solved at rest and the initial stresses are computed in each element. Then, for the portions to be excavated, the stresses in these elements are integrated in order to determine the nodal point forces which are applied in the opposite direction to the nodal points of ground–interface–liner system. The next step is to apply nodal forces to account for the weight of the liner. Unless otherwise stated in this section, the problem considered is a 9.0 m diameter tunnel with 40 cm concrete lining located at 18.0 m depth from the ground surface and the analyses are performed using both material and geometric nonlinearity and the interface simulation is taken into consideration. The water level is at the ground surface level and an undrained condition is assumed.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 5.
Finite element mesh of selected tunnel problem.
713
The soil was modeled using only the Drucker–Prager ultimate failure surface. This was done by fixing the initial cap position well beyond the maximum stress level that occurs in the problem and preventing the cap from contraction. The tunnel opening was inserted at different tunnel depths 2D, 3D, 4D, 5D noting that the tunnel depth is measured from the tunnel crown to the ground surface level. The analyses were performed considering two cases: the first case is the material nonlinearity only (MNO) while the second case is both material and geometric nonlinearity (MGN). It must be emphasized that the geometric nonlinearity is assumed throughout this parametric study taking into consideration large displacement, large rotations, but small strains. Fig. 7 shows the effect of the type of nonlinear analysis on both horizontal and vertical displacement at the ground surface level for different cases of tunnel depths. It can be observed that the surface settlement troughs approximate an inverted normal Gaussian distribution curve with a maximum value at the tunnel axis. On the other hand, the horizontal displacements have a zero value at the tunnel center line due to the symmetry of the problem and start increasing until they reach a maximum value near the point of inflection of the vertical settlement troughs. Also, it can be shown from the figure that the vertical and horizontal displacements decrease with the increase of tunnel depth. When performing a geometric nonlinearity analysis, both the vertical and horizontal displacements increased. At tunnel depth ⫽ 2D, the maximum vertical displacement increased by 4% whereas the maximum horizontal displacement increased by 3%. This indicates that the geometric nonlinearity analysis has a softening effect. A last conclusion which can be drawn from the figure is that the ratio of the maximum horizontal to vertical displacement ranges from 0.48 at tunnel depth ⫽ 2D to 0.39 at tunnel depth ⫽ 5D. 6.2. Effect of tunneling on the initial stresses of surrounding soil
Fig. 6.
Finite element mesh at contact surface.
6. Results of the parametric study 6.1. Effect of cavity only To study the effect of cavity only, a 5.0 m tunnel diameter opening was embedded in sandy soil which has a friction angle of 35° and cohesion (c) ⫽ 100.0 kN/m2.
The effect of the tunnel construction on the initial stress distribution in the surrounding soil mass is represented by the contour maps shown in Fig. 8. A negative sign indicates compressive stresses and the dimensions are plotted in meters. Fig. 8(a) and (b) show the vertical and horizontal normal stresses. It can be seen that the stress contours are affected all around the tunnel for a distance equal to its diameter in the case of vertical stresses. But in the case of horizontal stresses, the disturbed region increases to twice the diameter below the tunnel invert as shown in Fig. 8(b). Fig. 8(c) shows the contour map for shear stresses. It
714
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 7. Effect of type of nonlinear analysis on horizontal and vertical displacements at ground surface level for different cases of tunnel depths (cavity only).
should be noted that the stress values are small compared to the horizontal and vertical stresses and have bigger values near the tunnel. Stress values decrease radially with distance from the tunnel. Finally, Fig. 8(d) shows the contour map for maximum stresses, it is very similar to the contour map for horizontal stresses due to the fact that shear stresses are very small. 6.3. Effect of type of nonlinear analysis The effect of type of nonlinear analysis on the vertical and horizontal displacements at the ground surface level
is shown in Fig. 9. Both the horizontal and vertical displacements increase when using the material and geometric nonlinearity analysis. At tunnel depth ⫽ 2D, the maximum vertical displacement increased by 10% and the maximum horizontal displacement increased by 18%. At tunnel depth ⫽ 4D, the maximum vertical displacement increased by 9% and the maximum horizontal displacement increased by 22%. To illustrate the effect of the type of nonlinear analysis on the stress distribution in the tunnel lining, the normal stress distribution is plotted at the tunnel crown, springline, and invert as shown in Fig. 10. At the tunnel crown the case of MGN nearly did not alter the stress distri-
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
715
Fig. 8. Contour maps of stresses.
bution (only an increase of 0.5% in the maximum compressive stress). But at the springline, the maximum normal compressive stress was increased by 3.5% and at the tunnel invert by 5.7%. From Fig. 10, it can be noted that except for the assumption of linear elastic material the maximum tensile stress is very low due to the fact that concrete has
a low tensile strength compared to its compressive strength. A second remark which can be concluded from these figures is that after the region of tensile stresses reach their ultimate values, the normal stress distribution is released and becomes nearly zero, this represents zones of tensile cracks at the integration points in these sections.
716
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 9. Effect of nonlinear analysis on horizontal and vertical displacements at ground surface level for different cases of tunnel depths.
6.4. Effect of interface simulation Experimental results showed that slippage between concrete and soil is mostly dependent on the normal pressure at the contact surface and soil properties, i.e., its internal friction angle and cohesion (e.g., Clough and Duncan [7]). Based on these results, the static coefficient of friction, s, (see Section 3) can be taken as a ratio ranging from (0.65 tan–0.85 tan) where is the angle of internal friction of soil particles (since Boston Blue Clay has no cohesion). Therefore, s was assigned a value of 0.35. The effect of interface simulation on the ground sur-
face subsidence is shown in Fig. 11. The compatibility case is plotted in accordance with the interface one. They are nearly equal except in the case of tunnel depth ⫽ 2D, the interface condition slightly increases the magnitude of vertical displacements. The effect of interface simulation on the normal stress distribution at different tunnel sections can be shown from part (iv) of Fig. 10. It can be observed that taking interface condition into consideration has released some stresses in all the tunnel sections compared to the case of MGN analysis. While this reduction in stresses was 3.6% and 6.5% at the tunnel crown and tunnel springline, respectively, it was considerably smaller at the tunnel invert (about 1.0%).
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
717
Fig. 10. Normal stress distribution at: (a) crown section; (b) springline section; (c) invert section. (i) Case of linear elastic material; (ii) case of nonlinear elastic concrete model (MNO); (iii) case of nonlinear elastic concrete model (MGN); (iv) case of nonlinear elastic concrete model and interface simulation.
718
Fig. 11.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Effect of interface simulation on horizontal and vertical displacements at ground surface level for different cases of tunnel depths.
6.5. Effect of lining thickness
6.6. Effect of tunnel depth
To illustrate the effect of lining thickness, three different concrete lining thicknesses were considered (40, 50, 60 cm). The effect of changing the lining thickness on the ground subsidence is shown in Fig. 12. As expected, increasing the lining thickness results in a decrease in both the vertical and horizontal displacements. A similar observation can be made from Fig. 13 for the vertical displacements at the tunnel crown and inverts in addition to the horizontal displacements at the tunnel springline.
The effect of tunnel depth on the horizontal and vertical displacements at the ground surface level is shown in Fig. 9. Both figures show that increasing the tunnel depth results in a decrease in the amount of vertical displacements at the ground surface. Moreover, Fig. 9(b) shows that increasing the tunnel depth results in an increase in the width of the settlement trough. Fig. 13 illustrates the effect of tunnel depth on the inward vertical displacement at tunnel crown and the
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
719
Fig. 12. Effect of lining thickness on horizontal and vertical displacements at ground surface level for the case of tunnel depth ⫽ 2D.
tunnel invert, and the outward horizontal displacement at the springline. The results of the lining deformation indicate an elliptical mode of deformation. The figure shows that while increasing the tunnel depth results in an increase in the vertical displacement at the crown, it is associated with a decrease in the vertical displacement at the tunnel invert. As for the horizontal displacement at the springline, it increases with the increase in tunnel depth for the case of lining thickness ⫽ 40 cm while for lining thicknesses ⫽ 50 cm and 60 cm, the tunnel depth does not have a significant effect.
6.7. Effect of surcharge load
Fig. 14 shows the effect of surcharge load located in a width equals to half the tunnel diameter on the vertical displacement of ground surface at tunnel center line axis. The increase of vertical displacement with the increase of surcharge load is almost linear. MNO analysis and MGN analysis are nearly equal for smaller values of sur-
720
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 13.
Fig. 14.
Effect of tunnel depth on vertical displacement at crown, springline and invert for different cases of lining thickness.
Effect of surcharge load on vertical displacement at ground surface at tunnel center line for the case of tunnel depth ⫽ 2D.
charge load; afterwards, MGN results start to deviate from those of MNO to give larger values of displacements.
placements than the linear beam element. This follows from the fact that cracks in concrete result in a reduction in the stiffness.
6.8. Effect of lining modeling
7. Summary and conclusions
A comparison between the uncracked linear beam element and the nonlinear concrete quadrilateral planestrain element is given in Fig. 15. The beam element was modeled in the conventional way, i.e. its area and moment of inertia per meter length is input. Fig. 15 shows that the cracked concrete model gives larger dis-
The application of nonlinear soil plasticity and cracked concrete models together with finite element techniques to plane-strain tunnel problems taking into consideration the nonlinear nature of the contact between the concrete lining and the surrounding soil has been demonstrated.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
721
Fig. 15. Influence of lining modeling on horizontal and vertical displacements at ground surface for the case of tunnel depth ⫽ 2D.
Based on the numerical results outlined in Section 6, the following conclusions can be drawn: 1. The type of nonlinear analysis greatly affects the overall deformation pattern of the problem, either at the ground surface level or in the concrete lining. The material and geometric nonlinear analysis (MGN) has a softening effect on the deformations when assigned to soil elements. But when applied to concrete elements it has an insignificant effect since concrete is a relatively rigid material. 2. The normal stress distribution in the soil mass around the tunnel is affected to a distance approximately equal to the tunnel diameter. The shear stresses are
concentrated around the tunnel but they have very small values. 3. The vertical displacement profile at the ground surface level approximates an inverted normal gaussian distribution curve with a maximum value at the tunnel center line. The horizontal displacements at the ground surface reach their maximum value at the point of inflection of the vertical displacement profile. The ratio of the maximum horizontal to vertical displacement is about 0.39 to 0.48 depending on the tunnel depth. 4. The interface simulation does not alter the deformations at the ground surface significantly but it results in a release of the normal stresses in the tunnel
722
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
lining. Slippage usually occurs at angles making 45° with the tunnel axis. Tension release has not occurred anywhere around the tunnel circumference. 5. Increasing the lining stiffness leads to smaller deformations in the tunnel lining which leads, consequently, to much less displacement at the ground surface. 6. Deeper tunnel depths result in a decrease in the displacements at the ground surface and much wider settlement troughs. The inward vertical displacement at the tunnel crown and outward horizontal displacement at the tunnel springline have increased with increasing tunnel depth. However, the inward vertical displacement at the tunnel invert has decreased. 7. Modeling the tunnel lining by an uncracked linear beam element leads to less deformation. This is due to the fact that cracks in concrete result in a reduction of its stiffness. References [1] Bathe KJ. ADINA—A Finite Element Program for Automatic Dynamic Incremental Nonlinear Analysis. Watertown (MA): ADINA Engineering Inc., 1984.
[2] Bathe KJ. ADINA User’s Manual. Watertown (MA): ADINA Engineering Inc., 1984. [3] Bathe KJ. Finite Element Procedures in Engineering Analysis. New Jersey: Prentice-Hall, Inc., 1982. [4] Bathe KJ, Chaudhary A. A solution method for planar and axisymmetric contact problems. International Journal for Numerical Methods in Engineering 1985;21:65–88. [5] Bathe KJ, Ramaswamy S. On three-dimensional nonlinear analysis of concrete structures. Journal of Nuclear Engineering and Design 1979;52:385–409. [6] Chen WF, Baladi GY. Soil Plasticity—Theory and Implementation. Amsterdam: Elsevier, 1985. [7] Clough GW, Duncan JM. Finite element analysis of retaining wall behavior. Journal of the Soil Mechanics and Foundation Division, ASCE 1971;97(SM12):1657–73. [8] DiMaggio FL, Sandler IS. Material model for granular soils. Journal of the Engineering Mechanics Division, ASCE 1971;97(EM3):935–50. [9] Huang TK, Chen WF. Simple procedure for determining capplasticity-model parameters. Journal of Geotechnical Engineering, ASCE 1990;116(3):492–513. [10] Humphrey DN, Holtz RD. Cap Parameters for Clayey Soils. Innsbruck: Geomechanics, 1988:441-5. [11] Kupfer H, Hilsdorf HK, Rusch H. Behavior of concrete under biaxial stresses. Journal of the American Concrete Institute 1969;66(8):656–66. [12] McCarron WO, Chen WF. NFAP—User’s Manual, Structural Engineering Report No. CE-STR-86-4, School of Civil Engineering, Purdue University, Indiana, 1986.
Nonlinear analysis of tunnels in clayey/sandy soil with a concrete lining S.A.S. Youakim a, S.E.E. El-Metewally b, W.F. Chen b
c,*
a Arab Consulting Engineers, Cairo, Egypt Structural Engineering Department, El-Mansoura University, El-Mansoura, Egypt c School of Civil Engineering, Purdue University, West Lafayette, IN 47907, USA
Received 22 January 1998; received in revised form 12 January 1999; accepted 12 January 1999
Abstract In the nonlinear analysis of tunnels embedded in clayey or sandy soil with concrete lining, nonlinearity arises due to three parameters. The first parameter is the nonlinear material behavior of the lining and the surrounding soil. The second parameter is a result of the large deformations effect. The third parameter is the nonlinear nature of the contact between the concrete lining and the surrounding soil (essentially the slippage between the two materials). In this paper the three parameters of nonlinearity are modeled and implemented in a finite element program. 2000 Elsevier Science Ltd. All rights reserved. Keywords: Contact problems; Large displacements; Soil plasticity; Soil–structure interaction; Tunnels; Cracks
1. Introduction The response of an engineering problem such as soil– structure interaction subjected to static and/or dynamic loading is influenced significantly by the existence and behavior of the contact between the two different materials. Over the past years, analysis of soil–structure interaction problems has been carried out by assuming a complete bond at the interface, i.e., no relative motions are allowed. Although this assumption simplifies the analysis and design, it can introduce errors in the prediction of stresses and deformations. Moreover, in most of the latter researches in this field, a small deformation analysis has been assumed. But it was found that this assumption is not often true especially in the case of weak materials such as clayey soils. Another important aspect in the study of interaction problems is the nonlinear behavior of the foundation material and the surrounding soil mass. Although the application of classical limit analysis techniques to geotechnical problems gives rather satisfactory results, in many problems involving complex geometry or conditions in which knowledge of
* Corresponding author. Tel.: ⫹ 1-765-494-2254; fax: ⫹ 1-765496-1105.
deformations are important, the implementation of the advanced plasticity models may be very important. As a result of the rapid development of digital electronic computers with powerful numerical methods such as the finite element method, it is now possible to implement all the previously mentioned factors of nonlinearity in computer programs and perform a nonlinear finite element analysis for any selected problems and obtain results with a convenient degree of accuracy. The objective of this paper is to implement the nonlinear constitutive modeling of soil and concrete in addition to the assumption of large deformation analysis to a selected circular tunnel embedded in soft clayey soil or sandy soil and study how these factors may alter the predicted surface settlement or may change the stress distribution around the tunnel opening. Another intention of this paper is to study the phenomena of soil–structure interaction and how this may release or intensify the stresses in the tunnel lining. In this study, the finite element method is applied to a circular tunnel problem excavated by the shield method far away from the tunnel inlet and outlet faces. So, it is sufficient to consider the case of plane strain neglecting the advance velocity of the shield machine as it moves along the tunnel axis.
0141-0296/00/$ - see front matter 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 1 - 0 2 9 6 ( 9 9 ) 0 0 0 0 8 - 5
708
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
2. Nonlinear constitutive relations 2.1. Material modeling There exists a large variety of models which have been proposed in recent years to characterize the stress– strain and failure behavior of concrete and soil. All these models have certain inherent advantages and limitations which depend to a large degree on their particular application. In this paper the concrete material is modeled by a hypoelastic model whereas the soil is modeled by the cap model. It must be noted that in this paper our discussion is only limited to the time and temperature independent stress–strain constitutive relations. 2.1.1. Hypoelastic failure model for concrete Since concrete is an irreversible and load-path dependent material, it is appropriate to use the incremental stress–strain formulation, known as the hypoelastic model, in which the loading history effect is incorporated; thus, the load path-dependency behavior including the unloading behavior is represented. Bathe and Ramaswamy [5] developed a hypoelastic model for concrete based on a uniaxial stress–strain relation that is generalized to take biaxial and triaxial stress conditions into account. Tensile cracking and compression crushing conditions are identified using failure surfaces. The use of tensile and compression fail-
Fig. 1.
ure criteria including strain softening conditions prevents the unrealistically large stress and strain conditions predicted as in the case when using some plasticity models. With regard to compression failure, the failure envelope shown in Fig. 1 can be used. The shape of this failure surface is largely based on the experimental results by Kupfer et al. [11] after being extended to the three dimensional case. The compression failure envelope requires 24 discrete stress values to be input. Firstly, the values ip1/˜ c are input. These values define at what stress magnitudes tp1, the discrete two-dimensional failure envelopes for additional stress tp2 and tp3 are input. Then, the failure envelopes are defined by the failure stress value ijp3/˜ c (i ⫽ 1, . . ., 6; j ⫽ 1, 2, 3) that correspond to the stress magnitude tp2 ⫽ tp1, tp2 ⫽ tp3 ( is a constant), and tp2 ⫽ tp3, where tp1, t p2 and tp3 are the principal stresses at time t, with t p3 ⱕ tp2 ⱕ tp1. Tensile failure occurs if the tensile stress in a principal stress direction exceeds the tensile failure stress as expressed in Fig. 2. In this case, it is assumed that a plane of failure develops perpendicular to the principal stress direction. The effect of this material failure is that the normal and shear stiffnesses across the plane of failure are reduced, and the corresponding normal stress is released. Interested readers may refer to [5] for detailed derivation and development of the proposed model.
Triaxial compression failure envelope model.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
709
Fig. 2. Triaxial tensile failure envelope model.
2.1.2. Cap model for clayey soils The isotropic work-hardening plasticity cap model was first suggested by Drucker et al. cited in [9]. They introduced a spherical end cap to the Drucker–Prager model, in order to control the dilatancy of soil. DiMaggio and Sandler [8] proposed the concept of a strain hardening elliptical cap fitted to the perfectly-plastic (failure) surface of Drucker–Prager. An associated flow rule was employed and the movement of the cap is controlled by the increase or decrease of the plastic volumetric strain. From a theoretical point of view, the cap model is particularly appropriate to soil behavior, because it is capable of treating the conditions of stress history, stress
Fig. 3.
path dependency, dilatancy, and the effect of the intermediate principal stress. The schematic shape of this model is shown in Fig. 3. The loading function consists of two parts: 1. An ultimate failure envelope of the Drucker–Prager type which serves to limit the maximum shear stresses in the material and has the form: f 1 ⫽ ␣¯I1 ⫹ J1/2 2 ⫺ k ⫽ 0
(1)
where ¯I1 is the first invariant of the effective stress tensor; J2 is the second invariant of the deviatoric stress tensor; and ␣, k, are material constants related to the friction and cohesion of the soil.
Cap model; (a) I¯1 ⫺ J1/2 plane; (b) deviatoric plane. 2
710
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
2. A strain-hardening elliptic cap which has the form: f 2 ⫽ (I¯1 ⫺ L)2 ⫹ R2J2 ⫺ (X ⫺ L)2 ⫽ 0
(2)
where R is the ratio of the major to minor axis of the elliptic cap; and X and L define the ¯I1 value at the intersection of the elliptic cap with the ¯I1 axis and the failure function, respectively. The cap position, X, is controlled using the hardening function by DiMaggio and Sandler [8] as follows: X⫽⫺
冉 冊
1 ⑀pv ln 1 ⫺ D W
(3)
where W defines the maximum plastic volumetric compaction that the material can experience under hydrostatic loading (Fig. 4); D is a curve fitting parameter; and ⑀pv is the plastic volumetric strain. Detailed description of the cap model and the input parameters may be found elsewhere as in [6,9,10,12]. 2.2. Large deformation analysis In the present paper, the analyses are performed using Total-Lagrangian Formulation assuming both material and geometric nonlinearity considering large displacements, large rotations, but small strains. Bathe [3] gives detailed derivations of the governing finite element equations, stress–strain relations and the corresponding strain–displacement transformation matrices for different types of nonlinear analyses.
Fig. 4.
3. Constraint equations approach for soil–structure interaction
The response of a soil–structure interaction problem involving contacts or interfaces and subjected to static and/or dynamic loading is influenced significantly by the existence and behavior of the contacts. Many analyses of structure–foundation systems are carried out by assuming complete bonding at the interface; i.e., there is no relative motion. Although this assumption simplifies the analysis and design, it can introduce errors in the prediction of stresses and deformations. The main reasons for such errors are based on the fact that soil is neither homogeneous nor elastic and invariably much weaker than the structure it has to support. Moreover, experiments have shown that sliding between concrete/soil and steel/soil occurs at levels of interface shear which are significantly less than the shear limit of soil. Thus, analyses which assume a perfectly bonded interface condition over predict the shear transfer and depending on specific application, either over or under estimate structural response. Various methods have been proposed in the past in order to model the discontinuous behavior at the interface. Some workers have suggested the connection of the elements of soil and concrete to each other by discrete springs. Other applications have treated the interface by using elements of zero thickness or by thin-layered elements of small thickness. Some other methods use the approach of constraint equations or the total potential. This latter approach is adopted in this study and its main features are outlined in the following subsections.
Relationship between mean normal stress and volumetric strain.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
3.1. Concept
711
this case (tt is applied to equal dtn).
This solution approach is applicable to the analysis of contact problems involving sticking, frictional sliding and separation under large deformations. It was first presented by Bathe and Chaudhary [4]. It contains the following major aspects: 1. The total potential of the contact forces is included in the variational formulation to enforce the geometric compatibilities along the contact surfaces. 2. In the region of contact, surface tractions are evaluated from the externally applied forces. 3. The surface tractions between nodal points are employed to decide whether a nodal point is sticking, a sliding contact or releasing. 4. The number of equations due to the contact conditions is dynamically adapted to solve two equations for each node in contact if the node is in the sticking condition, and one equation if the node is in the sliding contact condition. This approach has the ability to solve contact problems whilst at the same time including large deformation effects, repeated contact and separation between the bodies at the interface. Hence, it will be considered in the present study. The approach considers two generic bodies denoted as contactor and target. In the finite element solution, the contactor contains the finite element boundary nodes that come into contact with the target segments or nodes. The basic conditions of contact along the contact surfaces are that no materials overlap can occur, and as a result, contact forces are developed and act along the regions of contact between the target and contactor. These forces are equal and opposite. The normal tractions can only excert compressive action, and the tangential tractions satisfy Mohr Coulomb’s law of friction. 3.2. Contact criteria Let tt ⫽ the developed tangential tractions along contact surfaces; tn ⫽ the compressive normal tractions; s ⫽ static coefficient of friction; and d ⫽ dynamic coefficient of friction (d ⱕ s). sticking contact:
The two bodies are assumed to be in sticking mode if 兩tt兩 ⱕ stn, which means that the frictional resistance during contact is sufficient to prevent sliding.
sliding contact:
sliding occurs between target and contactor bodies when 兩tt兩 > stn; i.e., the frictional capacity during contact is smaller than the applied tangential traction and
Tension contact:
The two bodies are assumed to separate from each other when tn ⬍ 0, and in this case both normal and/or tangential tractions are set to be zero.
4. Computer program—ADINA-96 ADINA-96 (Automatic Dynamic Incremental Nonlinear Analysis) is a computer program for the static and dynamic displacement and stress analysis of solids, structures and fluid–structure systems. The program has been written in FORTRAN 77, and can be employed to perform linear or nonlinear analysis in either plane strain, plane stress, axi-symmetric and general threedimensional analysis. ADINA-96 is written in a modular form so that any additional material models can be implemented without major programming effort. The program ADINA is the main module of the ADINA system which consists currently of the program ADINA for displacement and stress analysis, ADINAT for analysis of heat transfer and field problems, ADINA-IN for preparation and display of the input data and ADINA-PLOT for the display of the calculated solution results. The ADINA program was first developed by Professor K.J. Bathe [1]. Since then, it was subjected to various improvements in order to widen its element library, material models, and solution capabilities. The current version of the ADINA-96 program has two additional improvements added to the ADINA-84 computer program [2]. The first improvement is the implementation of the elliptic cap model for soils [10,12]. The elliptic cap model for soils was implemented after some modifications to justify the general control variables in the element library of ADINA-84 and the general analysis procedures. The second additional feature in ADINA-96 is the preprocessor program AUTOGEN [2] to assist in the preparation of the input file to ADINA by generating the model coordinates, two-dimensional element data, boundary conditions, and applied loads which balance the initial stresses in soil. In nonlinear analysis, ADINA-96 performs an incremental load-displacement analysis. After each increment of load, the computed displacement field is modified using an iterative procedure until an equilibrium configuration is reached.
5. Application to plane strain tunnel problem Soil response to construction procedures of tunnels have always been an area of great interest to engineers
712
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
since care must be taken to insure that structures adjacent to, or above, tunnel construction are not damaged by excessive ground movements. During tunnel excavation, the state of stresses in the ground continuum are significantly altered specially in soft soils. So, an immediate support must be provided to the surrounding ground; otherwise, the tunnel will collapse. On the other hand, prediction of the behavior of the tunnel lining during construction is a complicated problem due to the fact that neither the lining nor the ground around a tunnel behaves in an elastic manner. Another important factor is the nonlinear nature of the contact between the concrete lining and the surrounding soil mass. In this section, a parametric study is made to a selected circular tunnel problem concerning the type of nonlinear analysis, tunnel depth, tunnel thickness, interface behavior, lining modeling, and surcharge loads. 5.1. Input parameters for material models 5.1.1. Soil cap-model parameters A normally consolidated soft clayey soil, called the Boston-Blue Clay, is chosen to represent the soil around the tunnel cavity and its properties are presented in [6]. Humphrey and Holtz [10] used this type of soil in a nonlinear finite element study of reinforced embankments. The values of the different parameters of the capmodel for this type of soil can be found in [10,12]. 5.1.2. Concrete model parameters The behavior of the concrete lining of the tunnel is modeled by the hypoelastic model of Bathe and Ramaswamy [5] as discussed in Section 2.1.1. Table 1 summarizes the material and control parameters of the concrete model while Table 2 illustrates the parameters used to define the triaxial compression failure envelope, where E˜0, , ˜ t, ˜ c, e˜c, ˜ u and e˜u are the initial Young’s modulus, Poisson’s ratio, tensile stress, ultimate compressive stress, crushing stress, strain corresponding to crushing stress, ultimate compressive stress and ultimate compressive strain, respectively, and the tilde (~) stands for uniaxial condition.
Table 1 Material and control parameters for the concrete model Parameter
Value
E˜0 ˜ t ˜ c e˜c ˜ u e˜u
28 ⫻ 105 t/m2 (28.0 GPa) 0.18 400.0 t/m2 (4.0 MPa) 4000.0 t/m2 (40.0 MPa) 0.002 3500.0 t/m2 (35.0 MPa) 0.003
Table 2 Triaxial compression failure curves Principal stress ratios Curve number
p1/˜ c p3 at p2 ⫽ p1 p3 at p2 ⫽ 0.75 p3 p3 at p2 ⫽ p3
1
2
3
4
5
6
0.0 1.0 1.3 1.25
0.25 1.4 1.5 1.45
0.5 1.7 2.0 1.95
0.75 2.2 2.3 2.25
1.0 2.5 2.7 2.65
1.2 2.8 3.2 3.15
5.2. Problem geometry and sequence of analysis The chosen circular tunnel for this study is of 9.0 m outer diameter. This diameter is constant throughout the parametric study. In the case of a tunnel cavity without lining, a smaller diameter of 5.0 m is chosen in order to guarantee the stability of the opening, otherwise very large deformations would occur which may lead to a complete failure of the problem. Symmetry is assumed about the vertical axis so that one half the problem is considered in the analysis. The extent of the finite element mesh from the axis of symmetry is taken to be six times the diameter. This was found adequate to represent the ground subsidence because the deformations beyond this distance become nearly negligible. In the present analysis, it is assumed that the tunnel is excavated by the shield method and is located far from both the tunnel inlets and outlets, so that the assumption of plane-strain problem will be adequate. The soil continuum is represented by two dimensional 4–5 node quadrilateral elements. The concrete lining is represented by two dimensional 6-node quadrilateral elements. The contact between the lining and the soil is represented along the outer diameter of the tunnel. The nodes belonging to the soil and lining are initially assigned the same coordinates. Fig. 5 shows the whole geometry of the problem, while Fig. 6 demonstrates the finite element mesh at contact surfaces. The sequence of the analysis is organized as follows: first the problem is solved at rest and the initial stresses are computed in each element. Then, for the portions to be excavated, the stresses in these elements are integrated in order to determine the nodal point forces which are applied in the opposite direction to the nodal points of ground–interface–liner system. The next step is to apply nodal forces to account for the weight of the liner. Unless otherwise stated in this section, the problem considered is a 9.0 m diameter tunnel with 40 cm concrete lining located at 18.0 m depth from the ground surface and the analyses are performed using both material and geometric nonlinearity and the interface simulation is taken into consideration. The water level is at the ground surface level and an undrained condition is assumed.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 5.
Finite element mesh of selected tunnel problem.
713
The soil was modeled using only the Drucker–Prager ultimate failure surface. This was done by fixing the initial cap position well beyond the maximum stress level that occurs in the problem and preventing the cap from contraction. The tunnel opening was inserted at different tunnel depths 2D, 3D, 4D, 5D noting that the tunnel depth is measured from the tunnel crown to the ground surface level. The analyses were performed considering two cases: the first case is the material nonlinearity only (MNO) while the second case is both material and geometric nonlinearity (MGN). It must be emphasized that the geometric nonlinearity is assumed throughout this parametric study taking into consideration large displacement, large rotations, but small strains. Fig. 7 shows the effect of the type of nonlinear analysis on both horizontal and vertical displacement at the ground surface level for different cases of tunnel depths. It can be observed that the surface settlement troughs approximate an inverted normal Gaussian distribution curve with a maximum value at the tunnel axis. On the other hand, the horizontal displacements have a zero value at the tunnel center line due to the symmetry of the problem and start increasing until they reach a maximum value near the point of inflection of the vertical settlement troughs. Also, it can be shown from the figure that the vertical and horizontal displacements decrease with the increase of tunnel depth. When performing a geometric nonlinearity analysis, both the vertical and horizontal displacements increased. At tunnel depth ⫽ 2D, the maximum vertical displacement increased by 4% whereas the maximum horizontal displacement increased by 3%. This indicates that the geometric nonlinearity analysis has a softening effect. A last conclusion which can be drawn from the figure is that the ratio of the maximum horizontal to vertical displacement ranges from 0.48 at tunnel depth ⫽ 2D to 0.39 at tunnel depth ⫽ 5D. 6.2. Effect of tunneling on the initial stresses of surrounding soil
Fig. 6.
Finite element mesh at contact surface.
6. Results of the parametric study 6.1. Effect of cavity only To study the effect of cavity only, a 5.0 m tunnel diameter opening was embedded in sandy soil which has a friction angle of 35° and cohesion (c) ⫽ 100.0 kN/m2.
The effect of the tunnel construction on the initial stress distribution in the surrounding soil mass is represented by the contour maps shown in Fig. 8. A negative sign indicates compressive stresses and the dimensions are plotted in meters. Fig. 8(a) and (b) show the vertical and horizontal normal stresses. It can be seen that the stress contours are affected all around the tunnel for a distance equal to its diameter in the case of vertical stresses. But in the case of horizontal stresses, the disturbed region increases to twice the diameter below the tunnel invert as shown in Fig. 8(b). Fig. 8(c) shows the contour map for shear stresses. It
714
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 7. Effect of type of nonlinear analysis on horizontal and vertical displacements at ground surface level for different cases of tunnel depths (cavity only).
should be noted that the stress values are small compared to the horizontal and vertical stresses and have bigger values near the tunnel. Stress values decrease radially with distance from the tunnel. Finally, Fig. 8(d) shows the contour map for maximum stresses, it is very similar to the contour map for horizontal stresses due to the fact that shear stresses are very small. 6.3. Effect of type of nonlinear analysis The effect of type of nonlinear analysis on the vertical and horizontal displacements at the ground surface level
is shown in Fig. 9. Both the horizontal and vertical displacements increase when using the material and geometric nonlinearity analysis. At tunnel depth ⫽ 2D, the maximum vertical displacement increased by 10% and the maximum horizontal displacement increased by 18%. At tunnel depth ⫽ 4D, the maximum vertical displacement increased by 9% and the maximum horizontal displacement increased by 22%. To illustrate the effect of the type of nonlinear analysis on the stress distribution in the tunnel lining, the normal stress distribution is plotted at the tunnel crown, springline, and invert as shown in Fig. 10. At the tunnel crown the case of MGN nearly did not alter the stress distri-
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
715
Fig. 8. Contour maps of stresses.
bution (only an increase of 0.5% in the maximum compressive stress). But at the springline, the maximum normal compressive stress was increased by 3.5% and at the tunnel invert by 5.7%. From Fig. 10, it can be noted that except for the assumption of linear elastic material the maximum tensile stress is very low due to the fact that concrete has
a low tensile strength compared to its compressive strength. A second remark which can be concluded from these figures is that after the region of tensile stresses reach their ultimate values, the normal stress distribution is released and becomes nearly zero, this represents zones of tensile cracks at the integration points in these sections.
716
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 9. Effect of nonlinear analysis on horizontal and vertical displacements at ground surface level for different cases of tunnel depths.
6.4. Effect of interface simulation Experimental results showed that slippage between concrete and soil is mostly dependent on the normal pressure at the contact surface and soil properties, i.e., its internal friction angle and cohesion (e.g., Clough and Duncan [7]). Based on these results, the static coefficient of friction, s, (see Section 3) can be taken as a ratio ranging from (0.65 tan–0.85 tan) where is the angle of internal friction of soil particles (since Boston Blue Clay has no cohesion). Therefore, s was assigned a value of 0.35. The effect of interface simulation on the ground sur-
face subsidence is shown in Fig. 11. The compatibility case is plotted in accordance with the interface one. They are nearly equal except in the case of tunnel depth ⫽ 2D, the interface condition slightly increases the magnitude of vertical displacements. The effect of interface simulation on the normal stress distribution at different tunnel sections can be shown from part (iv) of Fig. 10. It can be observed that taking interface condition into consideration has released some stresses in all the tunnel sections compared to the case of MGN analysis. While this reduction in stresses was 3.6% and 6.5% at the tunnel crown and tunnel springline, respectively, it was considerably smaller at the tunnel invert (about 1.0%).
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
717
Fig. 10. Normal stress distribution at: (a) crown section; (b) springline section; (c) invert section. (i) Case of linear elastic material; (ii) case of nonlinear elastic concrete model (MNO); (iii) case of nonlinear elastic concrete model (MGN); (iv) case of nonlinear elastic concrete model and interface simulation.
718
Fig. 11.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Effect of interface simulation on horizontal and vertical displacements at ground surface level for different cases of tunnel depths.
6.5. Effect of lining thickness
6.6. Effect of tunnel depth
To illustrate the effect of lining thickness, three different concrete lining thicknesses were considered (40, 50, 60 cm). The effect of changing the lining thickness on the ground subsidence is shown in Fig. 12. As expected, increasing the lining thickness results in a decrease in both the vertical and horizontal displacements. A similar observation can be made from Fig. 13 for the vertical displacements at the tunnel crown and inverts in addition to the horizontal displacements at the tunnel springline.
The effect of tunnel depth on the horizontal and vertical displacements at the ground surface level is shown in Fig. 9. Both figures show that increasing the tunnel depth results in a decrease in the amount of vertical displacements at the ground surface. Moreover, Fig. 9(b) shows that increasing the tunnel depth results in an increase in the width of the settlement trough. Fig. 13 illustrates the effect of tunnel depth on the inward vertical displacement at tunnel crown and the
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
719
Fig. 12. Effect of lining thickness on horizontal and vertical displacements at ground surface level for the case of tunnel depth ⫽ 2D.
tunnel invert, and the outward horizontal displacement at the springline. The results of the lining deformation indicate an elliptical mode of deformation. The figure shows that while increasing the tunnel depth results in an increase in the vertical displacement at the crown, it is associated with a decrease in the vertical displacement at the tunnel invert. As for the horizontal displacement at the springline, it increases with the increase in tunnel depth for the case of lining thickness ⫽ 40 cm while for lining thicknesses ⫽ 50 cm and 60 cm, the tunnel depth does not have a significant effect.
6.7. Effect of surcharge load
Fig. 14 shows the effect of surcharge load located in a width equals to half the tunnel diameter on the vertical displacement of ground surface at tunnel center line axis. The increase of vertical displacement with the increase of surcharge load is almost linear. MNO analysis and MGN analysis are nearly equal for smaller values of sur-
720
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
Fig. 13.
Fig. 14.
Effect of tunnel depth on vertical displacement at crown, springline and invert for different cases of lining thickness.
Effect of surcharge load on vertical displacement at ground surface at tunnel center line for the case of tunnel depth ⫽ 2D.
charge load; afterwards, MGN results start to deviate from those of MNO to give larger values of displacements.
placements than the linear beam element. This follows from the fact that cracks in concrete result in a reduction in the stiffness.
6.8. Effect of lining modeling
7. Summary and conclusions
A comparison between the uncracked linear beam element and the nonlinear concrete quadrilateral planestrain element is given in Fig. 15. The beam element was modeled in the conventional way, i.e. its area and moment of inertia per meter length is input. Fig. 15 shows that the cracked concrete model gives larger dis-
The application of nonlinear soil plasticity and cracked concrete models together with finite element techniques to plane-strain tunnel problems taking into consideration the nonlinear nature of the contact between the concrete lining and the surrounding soil has been demonstrated.
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
721
Fig. 15. Influence of lining modeling on horizontal and vertical displacements at ground surface for the case of tunnel depth ⫽ 2D.
Based on the numerical results outlined in Section 6, the following conclusions can be drawn: 1. The type of nonlinear analysis greatly affects the overall deformation pattern of the problem, either at the ground surface level or in the concrete lining. The material and geometric nonlinear analysis (MGN) has a softening effect on the deformations when assigned to soil elements. But when applied to concrete elements it has an insignificant effect since concrete is a relatively rigid material. 2. The normal stress distribution in the soil mass around the tunnel is affected to a distance approximately equal to the tunnel diameter. The shear stresses are
concentrated around the tunnel but they have very small values. 3. The vertical displacement profile at the ground surface level approximates an inverted normal gaussian distribution curve with a maximum value at the tunnel center line. The horizontal displacements at the ground surface reach their maximum value at the point of inflection of the vertical displacement profile. The ratio of the maximum horizontal to vertical displacement is about 0.39 to 0.48 depending on the tunnel depth. 4. The interface simulation does not alter the deformations at the ground surface significantly but it results in a release of the normal stresses in the tunnel
722
S.A.S. Youakim et al. / Engineering Structures 22 (2000) 707–722
lining. Slippage usually occurs at angles making 45° with the tunnel axis. Tension release has not occurred anywhere around the tunnel circumference. 5. Increasing the lining stiffness leads to smaller deformations in the tunnel lining which leads, consequently, to much less displacement at the ground surface. 6. Deeper tunnel depths result in a decrease in the displacements at the ground surface and much wider settlement troughs. The inward vertical displacement at the tunnel crown and outward horizontal displacement at the tunnel springline have increased with increasing tunnel depth. However, the inward vertical displacement at the tunnel invert has decreased. 7. Modeling the tunnel lining by an uncracked linear beam element leads to less deformation. This is due to the fact that cracks in concrete result in a reduction of its stiffness. References [1] Bathe KJ. ADINA—A Finite Element Program for Automatic Dynamic Incremental Nonlinear Analysis. Watertown (MA): ADINA Engineering Inc., 1984.
[2] Bathe KJ. ADINA User’s Manual. Watertown (MA): ADINA Engineering Inc., 1984. [3] Bathe KJ. Finite Element Procedures in Engineering Analysis. New Jersey: Prentice-Hall, Inc., 1982. [4] Bathe KJ, Chaudhary A. A solution method for planar and axisymmetric contact problems. International Journal for Numerical Methods in Engineering 1985;21:65–88. [5] Bathe KJ, Ramaswamy S. On three-dimensional nonlinear analysis of concrete structures. Journal of Nuclear Engineering and Design 1979;52:385–409. [6] Chen WF, Baladi GY. Soil Plasticity—Theory and Implementation. Amsterdam: Elsevier, 1985. [7] Clough GW, Duncan JM. Finite element analysis of retaining wall behavior. Journal of the Soil Mechanics and Foundation Division, ASCE 1971;97(SM12):1657–73. [8] DiMaggio FL, Sandler IS. Material model for granular soils. Journal of the Engineering Mechanics Division, ASCE 1971;97(EM3):935–50. [9] Huang TK, Chen WF. Simple procedure for determining capplasticity-model parameters. Journal of Geotechnical Engineering, ASCE 1990;116(3):492–513. [10] Humphrey DN, Holtz RD. Cap Parameters for Clayey Soils. Innsbruck: Geomechanics, 1988:441-5. [11] Kupfer H, Hilsdorf HK, Rusch H. Behavior of concrete under biaxial stresses. Journal of the American Concrete Institute 1969;66(8):656–66. [12] McCarron WO, Chen WF. NFAP—User’s Manual, Structural Engineering Report No. CE-STR-86-4, School of Civil Engineering, Purdue University, Indiana, 1986.