Support pressure for circular tunnels in two layered undrained clay

Bibhsh Kumr, Jgdish Prsd Shoo

a Department of Civil Engineering, National Institute of Technology, Uttarakhand, 246174, India

b Department of Civil Engineering, Indian Institute of Technology, Roorkee, 247667, India

Keywords:Finite elements Layered clay Limit analysis Stability of tunnels

A B S T R A C T To estimate the required support pressure for stability of circular tunnels in two layered clay under undrained condition,numerical solutions are developed by performing finite element lower bound limit analysis in conjunction with second-order cone programming. The support system is assumed to offer uniform internal compressive pressure on its periphery.From the literature,it is known that the stability of tunnels depends on the overburden pressure acting over it,which is a function of undrained cohesion and unit weight of soil, and cover of soil. When a tunnel is constructed in layered undrained clay, the stability depends on the undrained shear strength,unit weight,and thickness of one layer relative to the other layer.In the present study,the solutions are presented in a form of dimensionless charts which can be used for design of tunnel support systems for different combinations of ratios of unit weight and undrained shear strength of upper layer to those of lower layer, thickness of both layers, and total soil cover depth.

1. Introduction

Tunnels and underground openings are constructed for better transport facilities, underground pipeline, canal and hydropower projects, etc. The stability of tunnels depends on the resistance of support system to prevent failure of soil mass. Thus, the support pressure required must be known prior to its design to maintain tunnel stability. In order to evaluate tunnel stability in homogeneous ground conditions, a number of investigations have been performed(Atkinson and Potts,1977;Mair,1979;Davis et al.,1980;Assadi and Sloan,1991;Wu and Lee,2003;Lee et al.,2006;Osman et al., 2006; Yang and Yang, 2010; Sahoo and Kumar, 2014a; Yang et al., 2015; Zhang et al., 2017). Atkinson and Potts (1977) have conducted centrifuge model tests assessing the stability of circular tunnels driven in sand supported by means of either compressed air or bentonite slurry. By conducting centrifuge model tests,Mair(1979) examined the response of tunnel driven in soft clay supported by compressed air pressure from inside of tunnel.Davis et al.(1980) computed the support pressure required for preventing collapse of underground openings/tunnels formed in soft clay under undrained condition based on the concept of limit analysis.Wu and Lee (2003) carried out centrifuge model tests to monitor ground movements and the associated collapse mechanism for both single and group of parallel tunnels in clay under undrained condition. Lee et al. (2006) conducted a series of centrifuge model tests and performed numerical analysis of these test results for evaluating the surface settlement trough,generation of excess pore water pressure and arching effect developed during the excavation process of single and parallel tunnels through soft clay.Osman et al.(2006) analyzed the stability of tunnel by examining the deformation pattern around unlined tunnel excavated in undrained clay using upper bound limit analysis. In order to determine the magnitude of required support pressure for tunnels located in various types of soils,solutions have been developed on the basis of finite element limit analysis with implementation of different mathematical programmings (Yang and Yang, 2010; Sahoo and Kumar, 2012, 2014a, 2018; Yang et al., 2015; Zhang et al., 2017).Limited studies are available in the literature to obtain the support pressure of tunnels located in the non-homogeneous soil (Sloan and Assadi, 1991; Hagiwara et al., 1999; Grant and Taylor, 2000;Nunes and Meguid, 2009; Wilson et al., 2011; Sahoo and Kumar,2013, 2019a). By performing finite element limit analysis and rigid block upper bound method, Sloan and Assadi (1991) and Wilson et al. (2011) computed the required support pressure for square and circular tunnels advanced in undrained cohesive soil taking into account linear variation of undrained strength with depth.For unlined circular tunnels in undrained clayey soil,Sahoo and Kumar(2013)performed stability analysis by employing finite element upper bound limit analysis considering a linear variation of undrained cohesion with depth. In layered soil, the ground movements induced by tunneling process depend on the stiffness and type of overlying soil mass. The effect of type and stiffness of overlying soil mass on the ground movements during tunneling was investigated by Hagiwara et al. (1999) with the help of centrifuge model tests for tunnels constructed in layered soil. By conducting a few centrifuge model tests, Grant and Taylor (2000)studied the stability of tunnels located in the purely cohesive soil underlain by granular soil. Nunes and Meguid (2009) carried out plane strain elastoplastic finite element analyses and laboratory tests to investigate the variation in bending stresses developed in lining of tunnel support in purely cohesive soil overlain by coarsegrained sand layer.Using lower bound finite element limit analysis in conjunction with second-order cone programming, Sahoo and Kumar (2019b) determined required support pressure to maintain the stability of tunnel advanced in undrained clay layer overlain by sand layer. Nevertheless, few researches seem to be available for predicting the required compressive pressure that should be offered by a support system of tunnel advanced in two layered clayey soil to maintain the tunnel stability.In the present study,it is aimed to produce solutions for estimating the compressive pressure required to support the soil surrounding the tunnel advanced in an undrained clay layer overlain by another undrained clay layer.The required support pressure is assumed to be acting normally and uniformly along the periphery of the tunnel. For the computations, finite element lower bound limit analysis and second-order cone programming were employed. The finite element limit analysis has been extensively used for solving various stability problems in geomechanics (Sahoo and Kumar, 2014b;Keawsawasvong and Ukritchon, 2017 a,b,c,d; Khuntia and Sahoo,2017). The influences of thickness of lower clay layer above tunnel crown,thickness of upper clay layer,strengths of both the layers in terms of undrained shear strength, and unit weights of both layers have been studied.The failure patterns,i.e.the proximity of stress states at failure, have also been presented for a few cases.

2. Problem definition and chosen domain

A circular tunnel having diameter D is excavated in undrained clay(φu1=0 )overlain by another layer of undrained clay(φu2=0 )at a cover depth of H below the horizontal ground surface, as illustrated in Fig.1.H1is the thickness of lower soil layer above the tunnel crown, and H2is the thickness of upper layer. To perform plane strain analysis, it is assumed that tunnel length is very large as compared to its diameter.Keeping the undrained shear strength/undrained cohesion (cu1) and unit weight (γ1) of the lower layer remain the same, the analysis was performed for the following cases: (i) undrained shear strength (cu1) and unit weight (γ1) of lower layer smaller than that of upper layer(cu2＞cu1and γ2＞γ1),and (ii) undrained shear strength and unit weight of lower layer greater than that of upper layer (cu2＜cu1and γ2＜γ1). The parameters c,γ and φ refer to the cohesion,unit weight and internal friction angle of soil mass, respectively, whereas the subscript u indicates the undrained condition, and 1 and 2 refer to the lower and upper layers,respectively.The undrained strengths of both the clay layers are modeled by employing Tresca yield criterion and an associated flow rule.In order to support soil mass surrounding the tunnel periphery for different cases as mentioned above, it is intended in the present study to calculate the required support pressure(σi)in the form of lining and associated anchorage system.The support pressure is expressed in the form of a dimensionless factor defined as σi/cu1.

Taking into account the symmetric nature of loading and geometry of the problem, only half of the domain as illustrated in Fig. 2a was used to perform the present analysis. The horizontal boundary (IJ) and vertical boundary (JK) were kept at a sufficient large distance from the centroid of tunnel by satisfying two conditions:(i)the plastic zone must be within the chosen domain;and(ii)the solutions will not change with further increase in the size of the domain.The location of the boundaries JK and IJ with Lh=3D-30D and Lv= 3D-20D has been found to be suffciient depending upon the values of H/D by performing sensitivity analysis as discussed in Section 5.

3. Analysis

The analysis was performed following the finite element formulation based on lower bound theorem of limit analysis developed by Sloan (1988) and the second-order conic optimization technique of Makrodimopoulos and Martin (2006). In the lower bound approach, the magnitude of the load at failure is determined in a statically admissible stress field, that is, the stress field satisfies the equilibrium and stress boundary conditions without violating yield criterion. For computing the magnitude of limit pressure required to maintain stability of the tunnel, the objective function is maximized subjected to a set of equality and inequality constraints imposed on the unknown nodal stresses.The equality constraints are produced in order to satisfy equilibrium of elements, statically admissible stress discontinuities, and stress boundary conditions; whereas, for the satisfaction of yield criterion, the inequality constraints are generated. The chosen domain as illustrated in Fig.2a is discretized into a number of three-noded triangular elements.None of the nodes are linked with more than a single element, and consequently, every interface between the adjacent elements becomes always the line of the stress discontinuity. It should be noted that each node remains unique to a particular element, and consequently, two or more nodes often share the same coordinates. There are three basic unknowns at each node, namely, horizontal normal stress (σxx), vertical normal stress(σyy),and shear stress(τxy).A typical finite element mesh for H/D = 6 with H1/D = 3 and H2/D = 3 is shown in Fig. 2c. All the meshes are generated in MATLAB (MathWorks, 2015) by writing programming codes.

For modeling the stress field under plane strain condition, a linear variation of stresses is chosen within each element and linear shape functions are employed. The stresses within each triangular element must satisfy the equilibrium conditions which are defined as follows:

Fig.1. Definition of problem.

Fig.2. (a)Chosen problem domain and associated stress boundary conditions;(b)Sign convention for stresses;and(c)A typical finite element mesh for H/D=6 with H1/D=3 and H2/D = 3. ‘+ve' shows that stresses are positive in these directions.

where γ is the body force per unit volume of soil mass(unit weight of soil mass) acting vertically in the downward direction. In Eqs.(1a)and(1b),normal tensile stresses are taken as positive and the positive direction of shear stresses and coordinate axes are shown in Fig.2b.Since each node is unique to an element,the unit weight of upper layer(γ2)and lower layer(γ1)are assigned at the nodes of elements of upper and lower layers, respectively, along the interface of two layers for satisfying the equilibrium condition given in Eq. (1b).

In order to obtain lower bound solutions close to the true solution, statically admissible stress discontinuities are included in the lower bound analysis (Drucker, 1953; Lysmer, 1970; Sloan,1988), and to obtain statically admissible stress discontinuity, the shear and normal stresses are required to be continuous along the common sides of any two adjacent elements.

No constraints on the stresses were imposed along the boundaries IJ and JK.Along boundary LK,the normal and shear stresses are made equal to zero while along the surfaces LM and NI,only shear stress is zero. Along the interfaces of the adjoining clayey soil and the tunnel lining, it has been specified that the developed shear stress is less than or equal to the undrained cohesion of the adjoining soil, that is,■■τxy

■■cu1.In order to obtain strictly lower bound solutions, the state of stresses should satisfy yield criterion everywhere within the problem domain. Under plane strain condition, the Tresca yield criterion for undrained clay may be expressed as

The inequality constraint (Eq. (2)) can be expressed as a set of second-order cone (Makrodimopoulos and Martin, 2006) by introducing a vector zi={z1iz2iz3i}Tat each node, satisfying the equation of second-order cone, i.e. z1iEq. (2) is represented by second-order cone constraints as

where σxx,i,σyy,iand τxy,iimply the horizontal normal stress,vertical normal stress and shear stress,respectively,at the i-th node;and N refers to the total number of nodes in the domain.

For modeling the support system, no exclusive elements were used. Following Sahoo and Kumar (2019b), the sides of triangular elements along the periphery of tunnel are treated as the segments of lining surface since for a given segment,the difference between chord length and arc length is negligible when the size of triangular elements is very small. The constraints developed for computing the uniform compressive support pressure to be exerted by the support system was already given in Sahoo and Kumar (2019b);however,for the sake of completeness,this part is again described herein. Fig. 3 shows a typical segment of the tunnel periphery,where σxx,1and σxx,2refer to the horizontal normal tractions in xdirection at nodes 1 and 2, respectively;σyy,1and σyy,2refer to the vertical normal tractions in y-direction at nodes 1 and 2, respectively; and τxy,1and τxy,2are the shear tractions on the vertical/horizontal planes at nodes 1 and 2,respectively.The angles θ1and θ2are the inclinations of radial lines passing through nodes 1 and 2,respectively, with the horizontal plane. The normal stresses σn,1and σn,2acting at nodes 1 and 2 on the tunnel periphery are expressed as

σn，1= σxx，1cos2θ1+σyy，1sin2θ1τxy，1sin(2θ1) (4a)

σn，2= σxx，2cos2θ2+σyy，2sin2θ2τxy，2sin(2θ2) (4b)

It is assumed that the internal support pressure (σi) acts uniformly on the tunnel periphery; hence, the following condition needs to be satisfied:

σi= σn，1= σn，2= … = σn，Nt(5)

where Ntis the total number of nodes along the tunnel periphery.

The lining used as support system is generally elastic in nature and considered as thin-walled tube under plane strain condition(Wood, 1975; Croll, 2001; Tamura and Hayashi, 2005; Wang and Koizumi, 2010; Wang et al., 2015). Though the required support pressure to make the surrounding soil mass in stable condition is not uniform along the tunnel periphery, to calculate the stiffness and thickness of tunnel lining based on the theory of thin-walled cylinders, the pressure distribution is assumed to be uniform,which can be calculated as follows:

where σ1and σ2are the longitudinal and hoop stresses in the lining,respectively;ε1and ε2are the longitudinal and hoop strains in the lining, respectively; and t, E and ν are the thickness, Young's modulus and Poisson's ratio of the lining, respectively.

For plane strain condition, ε1=0, thus we have

Lining stiffness per unit tunnel length is expressed as

The magnitude of limit support pressure at failure is optimized by satisfying the constraints imposed on the nodal stresses. For performing second-order cone programming in a statically admissible stress field, the objective function along with the constraints may be expressed as follows:

Objective function:

Maximize{σi} (15)

Constraints:

Assuming that tcis the change allowed in the thickness of lining of the cylinder, then we can obtain

Fig. 3. Various parameters associated with a very small length of segment of the tunnel periphery.Nodes 1 and 2 lie on the tunnel periphery and node 3 lies within the soil mass.The difference between chord length and arc length between nodes 1 and 2 is negligible as they are very close to each other.

where A is the global matrix of coefficients of all the constraints;Aequil, Abound, Adis, Asocpand Aucsare the coefficient matrices of constraints owing to element equilibrium, imposed boundary conditions, satisfying statically admissible stress discontinuities, a set of second-order cones and uniform compressive support pressure,respectively;I is the identity matrix;B is the global vector of constants;Bequil,Bbound,Bdis,Bsocpand Bucsare the global vectors of constants on account of element equilibrium, imposed boundary conditions, satisfying statically admissible stress discontinuities, a set of second-order cones and uniform compressive support pressure, respectively; σ is the global variable consisting of the nodal stresses, which can be defined as

z is the global auxiliary variable vectors defined as

It is worth mentioning that for determining the resistance to be offered by internal support pressure against active failure induced by overburden pressure and surcharge (if any) on the ground surface, the objective function is maximized. On the other hand, for computing the resistance to be provided by the overburden pressure and surcharge against the collapse caused by internal pressure which is termed as passive failure or blow out failure, the objective function is minimized, that is, Minimize{σi}.The pre-processing computer codes were written in MATLAB(MathWorks, 2015) and optimized solutions were developed on the basis of cone optimization using MOSEK (MOSEK ApS, 2015),an optimization toolbox available for MATLAB. To obtain the optimal lower bound solutions, the MOSEK optimizer is employed for solving different stability problems in geotechnical engineering (Khuntia and Sahoo, 2018; Sahoo and Khuntia, 2018;Ukritchon et al., 2018; Ukritchon and Keawsawasvong, 2018,2019).

4. Results and discussion

From the previous studies on homogeneous clayey soil(Assadi and Sloan, 1991; Osman et al., 2006; Wilson et al., 2011; Sahoo and Kumar, 2014a), it is noted that for a given diameter of the tunnel, the magnitude of required support pressure depends on the unit weight and cohesion of soil, and cover depth of the tunnel. In the present analysis, numerical solutions were obtained for a tunnel where soil mass lying over the tunnel is of a two-layered medium, that is, tunnel located in clay under undrained condition (φu1=0 ) with an overlay of relatively stiff or soft clayey soil under undrained condition (φu2= 0 ), where the unit weight and undrained cohesion of upper layer are greater or smaller than that of lower layer. The present solutions have been obtained for two cases:

(1) In the first case,the results are obtained for a homogeneous clayey soil with γ1D/cu1=γ2D/cu2equal to 1, 3 and 5;

(2) In the second case,keeping the unit weight and cohesion of the lower layer constant with γ1D/cu1equal to 1,3 and 5,the values of unit weight and cohesion of the upper layer are varied with:

(i) relatively higher cohesion(cu2＞cu1)for cu2/cu1=10,2.5 and 1.43 (or cu1/cu2= 0.1, 0.4 and 0.7) with relatively higher unit weight(γ2＞γ1)for γ2/γ1=1.2,1.4 and 1.6(or γ1/γ2= 0.83, 0.71 and 0.625), and

(ii) relatively lower cohesion(cu2＜cu1)for cu2/cu1=0.1,0.4 and 0.7 with relatively lower unit weight (γ2＜γ1) for γ2/γ1= 0.83, 0.71 and 0.625.

Fig.4. Variation of σi/cu1 with H2/D for H1/D = 1 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 1.

It is noted that the values of undrained cohesion and unit weight of soft to very stiff clay generally lie in the range of 14-22 kN/m3and 25-100 kPa, respectively. The computations were performed by choosing the values of(i)H/D ranging from 1 to 10,(ii)γ1D/cu1or γ2D/cu2varying from 1 to 5, and (iii) cu2/cu1between 0.1 and 10.

4.1. Variation of support pressure

In the present study, the support pressure is normalized with respect to the undrained cohesion of the lower layer(σi/cu1),which has been kept the same as that used for obtaining results in the case of homogeneous clay. In two-layered clay medium, the support pressure (σi) required for maintaining stability of tunnel depends on the thicknesses of upper layer (H2) and lower layer (H1) above the crown of tunnel, unit weight and cohesion of both layers. For various combinations of H1/D,cu2/cu1and γ2/γ1,the variation of σi/cu1with H2/D is shown in Figs. 4-12. It may be noted that the movement of soil mass towards tunnel at failure due to the action of gravitational load is resisted by the internal pressure offered by the support system. From the results presented in Figs. 4-12, the positive value of σiimplies that support system is required to be provided to prevent the collapse of tunnel; whereas, negative values of σiindicate that the tunnel remains stable without any support system.

Keeping the shear strength of lower layer constant, corresponding to particular values of H/D, H1/D, γ1D/cu1and γ2/γ1, the magnitude of required support pressure (i) decreases with increase in the undrained shear strength of upper layer when the upper layer soil is relatively stiffer compared to the lower one,that is, unit weight and shear strength of lower layer are smaller than that of upper layer(γ2＞γ1and cu2/cu1＞1)as shown in parts(a-c)in Figs. 4-12; and (ii) increases with reduction in the undrained shear strength of upper layer when the upper layer soil is relatively softer compared to the lower one, that is, the unit weight and shear strength of lower layer are greater than those of upper layer(γ2＜γ1and cu2/cu1＜1)as shown in parts(d-f)in Figs.4-12.It can be seen from Figs. 4-12 that the variation of σi/cu1is found to increase nonlinearly with increase in H2/D up to a certain limit,after which it increases linearly. It is known that the failure is mainly due to the weight of soil mass(overburden pressure)lying over the tunnel and the resistance to failure is offered by the undrained cohesion of soil mass. For smaller thickness of upper layer,the influence of driving force due to overburden pressure of upper soil layer may be less compared to the resistance owing to the cohesion of upper layer, leading to nonlinear variation of σi/cu1. On the other hand, with increase in the thickness of upper layer, the effect of driving force due to increase in overburden pressure of soil becomes dominant versus that of increase in the resistance owing to the cohesion of soil mass, which causes the stiffer or linear variation of σi/cu1.

Fig.5. Variation of σi/cu1 with H2/D for H1/D = 3 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 1.

Fig.6. Variation of σi/cu1 with H2/D for H1/D = 5 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 1.

Fig.7. Variation of σi/cu1 with H2/D for H1/D = 1 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 3.

Fig.8. Variation of σi/cu1 with H2/D for H1/D = 3 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 3.

Fig.9. Variation of σi/cu1 with H2/D for H1/D = 5 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 3.

Fig.10. Variation of σi/cu1 with H2/D for H1/D = 1 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 5.

For H2/D=5,the variations of σi/cu1with(i)cu2/cu1for different values of γ2/γ1with H1/D=5 and γ1D/cu1=3,(ii)γ2/γ1for different values of cu2/cu1with H1/D = 5 and γ1D/cu1= 3, (iii) γ1D/cu1for different values of γ2/γ1and cu2/cu1with H1/D=5,and(iv)H1/D for different values of γ2/γ1and cu2/cu1with γ1D/cu1=3 are shown in Figs. 13-16. It is observed that there is increase in the value of support pressure with(i)reduction in the undrained cohesion and increase in the unit weight of upper layer with respect to the lower layer for given values of H1/D,H2/D and γ1D/cu1,(ii)increase in γ1D/cu1for a given combination of H1/D,H2/D,γ2/γ1and cu2/cu1,and(iii)increase in H1/D for a given combination of H2/D, γ1D/cu1, γ2/γ1and cu2/cu1.

Fig.11. Variation of σi/cu1 with H2/D for H1/D = 3 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 = 0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 5.

Fig.12. Variation of σi/cu1 with H2/D for H1/D = 5 with different ratios of γ2/γ1: (a) cu2/cu1 =10, (b) cu2/cu1 = 2.5, (c) cu2/cu1 =1.43, (d) cu2/cu1 =0.7, (e) cu2/cu1 = 0.4, and(f) cu2/cu1 = 0.1, compared to homogeneous clay with γ1D/cu1 = γ2D/cu2 = 5.

From Figs. 4-12, it shows that compared to tunnel driven in homogeneous clay,the magnitude of σi/cu1required for the case of tunnel driven in two-layered clay medium(where the properties of lower layer remain the same as those of homogeneous clay) has been found to be dependent on the combined influence of H1/D,H2/D,γ1D/cu1, cu2/cu1and γ2/γ1, i.e.

(1) Greater only when cu2/cu1=1.43 for any values of H/D and γ1D/cu1, and dependent on the combined influence of H1/D,H2/D,γ1D/cu1,cu2/cu1and γ2/γ1for cu2/cu1equal to 2.5 and 10 for lower values of γ1D/cu1=1, as provided in Figs. 4-6.

(2) Smaller only when cu2/cu1= 0.7 for any values of H/D and γ1D/cu1, and dependent on the combined influence of H1/D,H2/D,γ1D/cu1,cu2/cu1and γ2/γ1for cu2/cu1equal to 0.4 and 0.1 for lower values of γ1D/cu1=1, as illustrated in Figs. 4-6.

(3) Always greater when cu2/cu1＞1 and always lower when cu2/cu1＜1 for higher values of γ1D/cu1=3 and 5, which can be noted from Figs. 7-12.

Further,for cu2/cu1＞1 in the case of layered soil,though the upper layer offers higher resistance than the homogeneous soil where the cohesion is the same as that of lower layer(cu1),the magnitude of σi/cu1corresponding to all the values of cu2/cu1＞1 for any combinations of H1/D, H2/D and γ1D/cu1is not always lower than that of homogeneous soil. This is due to the fact that the driving force causing failure is also higher for the tunnel driven in layered soil than that in homogeneous soil as γ2＞γ1.Similarly,when the resistance of upper layer is lower compared to that of homogeneous soil(cu2＜cu1),the magnitude of σi/cu1corresponding to all the values of cu2/cu1＜1 for any combinations of H1/D, H2/D and γ1D/cu1is not always greater than that of homogeneous soil because the driving force causing failure is also lower for the tunnel driven in layered soil in comparison to the tunnel driven in homogenous soil as γ2＜γ1.

Fig.13. Variation of σi/cu1 with cu2/cu1 for (a) cu2/cu1＞1 and γ2/γ1 ＞1 and (b) cu2/cu1＜1 and γ2/γ1 ＜1, for H1/D = 5, H2/D = 5 and γ1D/cu1 =3.

Fig.14. Variation of σi/cu1 with γ2/γ1 for (a) cu2/cu1＞1 and γ2/γ1 ＞1, and (b) cu2/cu1 ＜1 and γ2/γ1 ＜1, for H1/D = 5, H2/D = 5 and γ1D/cu1 =3.

Fig.15. Variation of σi/cu1 with γ1D/cu1 for (a) cu2/cu1 =10 and γ2/γ1 ＞1 and (b) cu2/cu1 =0.1 and γ2/γ1 ＜1, for H1/D = 5 and H2/D = 5.

4.2. Failure patterns

The proximity of the stress state to shear failure, that is, the failure pattern,has been examined from the stress state obtained at the element centroid and expressed in dimensionless term,namely,t/s, where t =and s = 2cu, where cubecomes equal to cu1and cu2for the elements in the lower and upper layers, respectively. A point within the chosen domain will be in a state of shear failure when the ratio t/s tends to unity. The failure patterns are generated in a manner such that the color of an element becomes darker when approaching the shear failure (t/s = 1). The failure patterns have been generated for H/D = 6 (H1/D=3 and H2/D=3)and γ1D/cu1=3,as shown in Fig.17 in the case of layered clay for different values of cu2/cu1and γ2/γ1keeping the cohesion of lower layer constant.

Fig.16. Variation of σi/cu1 with H1/D for (a) cu2/cu1 =10 and γ2/γ1 ＞1, and (b) cu2/cu1 =0.1 and γ2/γ1 ＜1, for H2/D = 5 and γ1D/cu1 =3.

Fig.17. Proximity of shear failure for layered clay with H1/D=3,H2/D=3 and γ1D/cu1=3:(a)cu2/cu1=10 and γ2/γ1=1.6;(b)cu2/cu1=1.43 and γ2/γ1=1.6;and(c)cu2/cu1=0.7 and γ2/γ1 =0.625.

Fig.18. Variation of σi/cu1 with Lh/D keeping Lv/D = 15 for H1/D = 3, H2/D = 3, γ1D/cu1 =3, cu2/cu1 =10 and different values of γ2/γ1.

From Fig. 17a, for a given value of γ2/γ1in the case of twolayered clayey soil, it can be noted that the yielding of soil mass depends upon the magnitude of undrained cohesion of both layers relative to each other, that is, the zone of yielding is larger as expected in the lower layer where the cohesion is relatively lower compared to that of upper layer (cu2/cu1= 10). By comparing Fig.17a,b,it is observed that the size of failure zone tends to reduce in the lower layer and increase in the upper layer with the reduction in the values of cohesion of upper layer. Further, though the cohesion of upper layer is less, the zone of yielding is found to be smaller in the upper layer for cu2/cu1= 0.7 as shown in Fig. 17c,when compared with that of cu2/cu1=1.43 as shown in Fig.17b. It needs to be mentioned that the zone of failure in the upper layer also depends on the ratio of unit weight of both layers, that is, the unit weight of upper layer is smaller for the case shown in Fig.17c than that presented in Fig.17b.

4.3. Remarks

For H1/D=3,H2/D=3,γ1D/cu1=3,cu2/cu1=10 and γ2/γ1=1.2,1.4 and 1.6,the variation of σi/cu1with Lh/D is presented in Fig.18.It could be noted that the magnitude of σi/cu1increases continuously with increase in domain size up to a certain limit (Lh/D is approximately equal to 12), beyond which with further increase in the domain size, the values of σi/cu1remains constant. Further, from Fig.17a,the developed plastic zone is observed to be well contained within the chosen domain with Lhand Lvequal to 25D and 15D,respectively.Thus,the chosen domain with Lh=25D and Lv=15D for determining the magnitude of σi/cu1when H1/D = 3, H2/D = 3,γ1D/cu1=3, cu2/cu1=10 and γ2/γ1=1.6 is proved to be sufficient.Similarly, for different combinations of H1/D, H2/D, γ1D/cu, cu2/cu1and γ2/γ1, the domain size was decided in the present analysis.

Fig.19. Comparison of present solutions with (a)numerical solutions of Wilson et al.(2011)and (b) experimental results of Mair (1979) for tunnel located in homogenous clayey soil.

5. Comparison of present solutions with those available in the literature

Since no solutions seem to be available in the literature for determining the required support pressure to maintain the tunnel stability in two-layered clayey soil,for the purpose of validation of present analysis,the variation of σi/cuwith H/D for different values of γD/cuobtained from the present analysis is compared with the numerical solutions of Wilson et al. (2011) and centrifuge experimental results of Mair(1979)for a tunnel located in homogeneous clay under undrained condition. The comparisons are provided in Fig.19a and b. In these figures, cuand γ are the undrained shear strength and unit weight of homogenous undrained clay(cu= cu1= cu2, and γ =γ1=γ2), respectively. Wilson et al. (2011)analyzed the stability of tunnel using finite element lower bound limit analysis in conjunction with nonlinear programming. The solutions obtained from these two analyses are found to be almost merging with each other. The present results also match closely with the experimental results of Mair(1979).

6. Conclusions

The support pressure (σi) required for the stability of a circular tunnel formed in clay layer overlain by relatively stiffer and softer clay layer has been computed by using lower bound finite element limit analysis in combination with the second-order conic programming.The analysis was performed considering the properties of both layers corresponding to undrained condition.The required support pressure is presented in terms of normalized normal compressive stress,defined as σi/cu1,and the variations of σi/cu1for various combinations of H1/D,H2/D,γ1D/cu1,cu2/cu1and γ2/γ1have been established in the form of dimensionless charts,which can be used for the purpose of design of tunnel support system in a twolayered undrained clayey soil. The following conclusions can be drawn from the present study:

(1) For a particular combination of H1/D,H2/D,γ1D/cu1and γ2/γ1,the magnitude of σi/cu1is found to be continuously reducing with increase in the undrained cohesion of the upper layer when tunnel is located in the lower layer of relatively smaller stiffness than that of the upper layer(cu2/cu1＞1 and γ2＞γ1),and continuously increasing with reduction in the undrained cohesion of the upper layer when tunnel is driven in the lower layer of relatively greater stiffness than that of the upper layer (cu2/cu1＜1 and γ2＜γ1).

(2) With increase in the unit weight of upper layer with respect to the lower layer for given values of H1/D,H2/D and γ1D/cu1,the magnitude of required support pressure increases.

(3) The increase in support pressure for a particular value of tunnel cover depth is nonlinear for lower values of H2/D and becomes linear for higher values of H2/D. The rate of nonlinearity is found to be higher when the undrained cohesion of upper layer is relatively greater than that of lower layer.Furthermore,the support pressure also increases with increase in the thickness of lower layer when the thickness of upper layer is constant.

(4) With increase in the values of γ1D/cu1, the magnitude of required support pressure also increases when the thicknesses of both layers, ratio of unit weight of upper to lower layer, ratio of cohesion of upper layer to lower layer remain constant.

(5) When the support pressure of a tunnel driven in layered soil with cu2＞cu1and γ2＞γ1is compared with that of a tunnel driven in homogeneous soil having the same properties as that of lower layer (cu1and γ1), it has been noted that the magnitude of support pressure required is not always lower for a tunnel driven in a layered soil,although the cohesion of upper layer is higher than that of homogeneous soil. This is because the driving force inducing failure is also higher for the tunnel driven in layered soil than that in homogeneous soil as γ2＞γ1.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.