Examples • Verification Problems
Lined Circular Tunnel in an Elastic Medium with Anisotropic Stresses
Problem Statement
Note
To view this project in 3DEC, use the menu command . Choose “BlockSel/ VerificationProblems/ linedCircularTunnel” and select “LinedCircularTunnel.prj” to load. The main data file used is shown at the end of this example. The remaining data files can be found in the project.
A circular tunnel with a radius of 5 m is located at 30 m depth in a soft elastic soil (\(E\) = 48 MPa, \(\nu\) = 0.34, \(\rho\) = 2000 kg/m^{3}). The insitu stresses are 600 kPa vertical and 300 kPa horizontal. The tunnel is supported by a 125 mm thick shotcrete liner (\(E=25 GPa\), \(\nu=0.15\)). We assume that the support is installed simultaneously with the excavation in a preexisting anisotropic biaxial stress field. The support displacements and internal stresses (in terms of axial thrust and moment), and the interface contact stresses are computed under planestrain conditions for the two limiting conditions of noslip (no relative shear displacement) and fullslip (no shear stress transmission) at the groundsupport interface. These computed values are compared with the analytical solution of Einstein and Schwartz (1979).
Analytical Solution
The analytical solution (Einstein and Schwartz 1979) is expressed using the notation in Figure 1. The support displacements consist of a radial, \(u_s\), and a tangential, \(v_s\), component. The internal stresses consist of an axial thrust, \(T\), and a bending moment, \(M\). The interface contact stresses consist of a normal, \(\sigma_R\), and a shear, \(\tau_{R\theta}\), component.
NoSlip Solution — For the noslip solution, the interface boundary condition consists of no relative shear displacement between the ground and the support. The noslip solution is given in the following equations:
where: 
\(\theta\) 
= 
angular location (counterclockwise with respect to horizontal); 
\(R\) 
= 
tunnel radius; 

\(P\) 
= 
vertical stress; 

\(K\) 
= 
ratio of horizontaltovertical stress; 

\(E\) 
= 
Young’s modulus of the ground mass; 

\(\nu\) 
= 
Poisson’s ratio of the ground mass; and 

\(a^*_0,a^*_2,b^*_2\) 
= 
dimensionless coefficients (see the next equations below). 
where: 
\(C^*\), \(F^*\) 
= 
compressibility and flexibility ratios, respectively (see the next two equations). 
where: 
\(E_s\) 
= 
Young’s modulus of the support; 
\(\nu_s\) 
= 
Poisson’s ratio of the support; 

\(A_s\) 
= 
average crosssectional area of support per unit length of tunnel 

(for support of constant thickness \(t\), \(A_s = t\)); and 

\(I_s\) 
= 
moment of inertia of support per unit length of tunnel 

(for support of constant thickness \(t\), \(I_s = t^3/12\)). 
FullSlip Solution — For the fullslip solution, the interface boundary condition consists of no shear stress transmission between the ground and the support. The fullslip solution is given in the following equations:
where: 
\(\theta\) 
= 
angular location (counterclockwise with respect to horizontal); 
\(R\) 
= 
tunnel radius; 

\(P\) 
= 
vertical stress; 

\(K\) 
= 
ratio of horizontaltovertical stress; 

\(E\) 
= 
Young’s modulus of the ground mass; 

\(\nu\) 
= 
Poisson’s ratio of the ground mass; and 

\(a^*_0\), \(a^*_2\) 
= 
dimensionless coefficients (see the next equations). 
where: 
\(C^*\), \(F^*\) 
= 
compressibility and flexibility ratios, respectively (see equation (1) above). 
3DEC Model
The problem is described in terms of the following parameters from the analytical solution:
Geometry 

tunnel radius (\(R\)) 
5 m 
InSitu Stresses 

vertical stress (\(P\)) 
600 kPa 
ratio of horizontaltovertical stress (\(K\)) 
0.5 
(\(\sigma_{xx}=300\) kPa; \(\sigma_{zz}=600\) kPa) 

Ground Mass Properties 

Young’s modulus (\(E\)) 
48 MPa 
Poisson’s ratio (\(\nu\)) 
0.34 
(shear modulus = 17.91 MPa; bulk modulus = 50 MPa) 

Linear Properties 

Young’s modulus (\(E_s\)) 
25 GPa 
Poisson’s ratio (\(\nu_s\)) 
0.15 
thickness (\(t\)) 
125 mm 
The 3DEC model simulates a thin slice of a circular tunnel in an infinite elastic ground mass with a preexisting anisotropic biaxial stress field subjected to planestrain conditions. The geometry of the 3DEC model is shown in Figure 2. The farfield boundaries are placed at a distance of 20 times the tunnel radius to approximate infinite boundaries. The insitu stresses are installed in all zones and also applied as loads acting on the farfield boundaries. Planestrain conditions are enforced by including a thin slice of material in the \(y\)direction and imposing symmetry boundary conditions on these two surfaces. Symmetry boundary conditions are also imposed on the planes at \(x=0\) and \(z=0\). For gridpoints, this requires maintaining zero displacement normal to the plane. For nodes, this requires maintaining zero displacement normal to the plane and zero rotation about two axes that lie in the plane. (The appropriate nodal conditions are imposed by realigning and then fixing the appropriate nodelocal systems, and also specifying the proper velocityfixity conditions — see Advancing Lined Tunnel for a description of the procedure.)
The 3DEC grid contains a single layer of zones in the \(y\)direction and is graded as one moves away from the tunnel (see Figure 2). The model resolution is 24 zones along the tunnel boundary. Liner structural elements with a crosshatch mesh pattern (see LinerReinforced Beam for reasons why this mesh pattern is chosen) are attached to the zone faces lying along the tunnel boundary (see Figure 3). The linerzone interface stiffnesses (\(k_n\) and \(k_s\)) are chosen using the linerstiffness equation, and increasing the value by a factor of 100 as suggested in the text following this equation. (We will confirm below that the criterion of small interface deformation is met for our system.)
Results and Discussion — Qualitative Assessment
We confirm that the interface deformation is small relative to the zone deformation (and thus confirm that the elastic interface stiffnesses, \(k_n\) and \(k_s\), are large enough) by plotting the displacements of the gridpoints and nodes at the tunnel crown and springline (see Figure 4). The relative displacement between the gridpoints and nodes is small compared with the gridpoint displacement.
If the model is run with no support (by deleting the liner elements), the tunnel crown and springline both move inward (see Figure 5). When the support is included, the tunnel crown still moves inward, but the tunnel springline moves outward (see Figure 6), because the support resists the inward ground movement. If the liner is allowed to slip at the supportground interface (by setting the liner property of couplingcohesionshear equal to zero), a relative shearing motion occurs between the liner and the ground (see Figure 7).
The axial thrust, \(T\), and bending moment, \(M\), in the liner for both the noslip and fullslip cases are shown in Figure 8 to Figure 11. When the liner is allowed to slip, the axial thrust becomes more uniform and the bending moment increases slightly. We define a liner surface coordinate system whose \(x\)axis lies along the tunnel axis (in the global \(y\)direction), and whose \(z\)axis is normal to the shell midsurface (pointing inward). This corresponds to the surface = (0, 1, 0) direction during the stress recovery procedure. (The surface coordinate system needs to be specified. See Stress Recovery Procedure for details.) In terms of this system, the axial thrust corresponds with the membrane stress resultant, \(N_y\), and the bending moment corresponds with the bending stress resultant, \(M_y\). In these plots, the liner is oriented such that lefttoright corresponds with the angular location \(\theta\) in the closedform solution varying from zero to ninety degrees. Note that these values of \(N_y\) and \(M_y\) are of opposite sign to the values of \(T\) and \(M\) from the analytical solution. (Positive \(N_y\) extends the shell, and positive \(M_y\) produces positive stress at the outer fiber on the side of the shell defined by the positive \(z\)direction of the surface system.)
The interface contact stresses, \(\sigma_R\) and \(\tau_{R\theta}\), for both the noslip and fullslip cases are shown in Figure 12 to Figure 15. These values are obtained from the normal and shear coupling springs that join the liner nodes to the zones. Positive normal stresses indicate separation, and this convention is opposite to that of the analytical solution. The plotted shear stresses depict magnitude only; the direction can be printed or accessed by the FISH functions struct.liner.normal.dir and struct.liner.shear.dir.
break
A quantitative comparison of the support displacements, the internal stresses, and the interface contact stresses with the analytical solutions for both the noslip and fullslip cases is provided in Figure 16 to Figure 24. The FISH functions extract the necessary values from the 3DEC model and compare them with the analytical values. Either the noslip or fullslip analytical solution is used by setting the no_slip FISH variable to 1 or 0, respectively. The support displacements (\(u_s\) and \(v_s\)) are sampled at the nodes using the struct.node.disp.global function. The internal stresses (\(T\) and \(M\)) are sampled at the element centroids using the struct.shell.resultant function. Also, the interface contact stresses (\(\sigma_R\) and \(\tau_{R\theta}\)) are sampled at the linerSEL centroids by computing the average value from the three coupling springs at the linerSEL nodes using the struct.liner.normal.stress and struct.liner.shear.stress functions.
For both the noslip and fullslip cases, the responses compare well with the analytical solution. The support displacements match the analytical solution, with a maximum discrepancy of approximately 2% at the tunnel crown. The axial thrusts match the analytical solution, with a maximum discrepancy of approximately 1.2%. The interface contact stresses match the analytical solution away from the modeled ends (at \(\theta\) equals 0 and 90 degrees). There is an error of approximately 10% at the end locations (note the outliers in Figure 19 and Figure 24) that arises from the lack of symmetry of the mesh. The area assigned to the two nodes at each end is not the same because it is the summation of onethird of the area of each element using the node. At each end, one node has twice the area of the other, and this affects the stress computation such that the result is nonsymmetric. Notice that the two outliers bound the analytical solution, and that their average value equals the analytical solution.
Reference
Einstein, H. H., and C. W. Schwartz. “Simplified Analysis for Tunnel Supports,” J. Geotech. Engr. Div., 105(GT4): 499518 (1979).
break
Data File
LinedCircularTunnel.dat
model new
model largestrain off
fish automaticcreate off
[global t1 = 'Lined Circular Tunnel ']
[global t2 = 'in an Elastic Medium with Anisotropic Stresses']
[global t = t1 + t2]
model title [t]
; Create the ground mass (zones).
block create tunnel blockradial 40 blockstangential 12 blocksaxial 1 ...
boundary 20 length 0 0.4 radius 5.0 radiusratio 1.05
block delete range posx 105 0
block delete range posz 105 0
block join
block zone generate edgelength 10.0
block group 'rock'
block group 'tunnel' range cylinder end1 0 0 0 end2 0 0.4 0 radius 5.0
block face group 'surface' range groupintersection 'rock' 'tunnel'
block delete range group 'tunnel'
; Material model and properties
block zone cmodel assign elastic
block zone property bulk 5e7 shear 1.791e7 dens 2500
; Create the support (linerSELs).
struct liner create byblockface range group 'surface'
struct liner cmodel assign elastic
struct liner property young 2.5e10 poisson 0.15 thickness 0.125 ...
couplingstiffnessnormal 2.0763e10 ...
couplingstiffnessshear 2.0763e10 ...
couplingcohesionshear 1e20
; Install insitu stresses in entire grid.
block zone initialize stress xx 3e5 yy 3e5 zz 6e5
; Specify boundary conditions.
; Apply stresses at farfield boundaries
block face apply stressxx 3e5 range posx 105.0
block face apply stresszz 6e5 range posz 105.0
; For the grid points (symmetry conditions):
block gridpoint apply velocityx 0.0 range posx 0.0
block gridpoint apply velocityz 0.0 range posz 0.0
block gridpoint apply velocityy 0.0 range posy 0.0
block gridpoint apply velocityy 0.0 range posy 0.4
; For the nodes (symmetry conditions):
struct node systemlocal x (1,0,0) y (0,1,0) ...
range positionx 0 ; x=0 plane
struct node fix systemlocal ...
range positionx 0
struct node fix velocityx rotationy rotationz ...
range positionx 0
struct node systemlocal x (0,0,1) y (0,1,0) ...
range positionz 0 ; z=0 plane
struct node fix systemlocal ...
range positionz 0
struct node fix velocityx rotationy rotationz ...
range positionz 0
struct node fix velocityy rotationx rotationz ...
range positiony 0 ; y=0 plane
struct node fix velocityy rotationx rotationz ...
range positiony 0.4 ; y=Ymax plane
; Histories
model history mechanical ratiolocal
block history displacementx position (5,0,0) ; tunnel springline
block history velocityx position (5,0,0)
struct node history displacementx position (5,0,0)
struct node history velocityx position (5,0,0)
block history displacementz position (0,0,5) ; tunnel crown
block history velocityz position (0,0,5)
struct node history displacementz position (0,0,5)
struct node history velocityz position (0,0,5)
history interval 10
;
model save 'starting'
; No support case
struct liner delete
model solve
model title [t]
model save 'No Support'
; Slip case
model restore 'starting'
struct liner property couplingcohesionshear 0
model solve
model title [t + "  Full Slip"]
model save 'Full Slip'
; No slip case
model restore 'starting'
model solve
model title [t + "  No Slip"]
model save 'No Slip'
program return
Was this helpful? ...  Itasca Software © 2023, Itasca  Updated: Sep 13, 2023 