# Cylindrical Hole in an Infinite Mohr-Coulomb Medium

Problem Statement

Note

To view this project in 3DEC, use the menu command Help ► Examples…. The main data files used are shown at the end of this example. All data files are found in the project.

This problem concerns determination of stresses and displacements around a cylindrical hole in an infinite elasto-plastic medium subjected to a compressive, isotropic far-field stress. The objective of this problem is to test the plasticity calculations in 3DEC. The medium is assumed to be linearly elastic, perfectly plastic. The Mohr-Coulomb criterion is used as the plastic limit; the material deforms plastically according to the associated flow rule (i.e., the dilation angle is equal to the friction angle).

The material is assigned the following properties:

density ($$\rho$$): 2600 kg/m3

shear modulus ($$G$$): 2.8 GPa

bulk modulus ($$K$$): 3.9 GPa

cohesion ($$c$$): 3.45 MPa

friction angle ($$\phi$$): 30°

dilation angle ($$\psi$$): 30°

An isotropic far-field stress is equal to 30 MPa. The radius of the hole is 1 m, and is assumed to be small compared to the length of the cylinder. (This permits the use of the plane-strain solution for comparison.)

Analytic Solution

The radius of the plastic zone, $$R_o$$, around a circular cylindrical hole in an infinite Mohr-Coulomb material (Figure 1) is given by the analytical solution derived by Salençon (1969):

$R_o = a\ \left ( \thinspace {2 \over {K_p + 1}}\ {P_o + {q \over {K_p - 1}} \over {P_i + {q \over {K_p - 1}}}}\thinspace \right )^{1/(K_p - 1)}$

where:

$$a$$ = radius of hole;

$$K_p$$ = $$(1 + \sin \phi)/ (1 -\sin \phi)$$;

$$q$$ = $$2c \tan (45 + \phi /2)$$;

$$P_o$$ = initial in-situ stress; and

$$P_i$$ = internal pressure.

According to the solution, the radial stress at the elastic/plastic interface is

$\sigma_{re} = {1 \over {K_p + 1}}\ (2 P_o - q)$

The stresses in the plastic zone are:

$\sigma_r = -{q \over {K_p - 1}} + \left(P_i + {q \over {K_p - 1}}\right) \cdot \left({r \over {a}}\right)^{K_p - 1}$
$\sigma_\theta = -{q \over {K_p - 1}} + K_p \left(P_i + {q \over {K_p - 1}}\right) \cdot \left({r \over {a}}\right)^{K_p - 1}$

where $$r$$ is the distance to the center of the hole; and the stresses in the elastic zone are:

$\sigma_r = P_o - (P_o - \sigma_{re})\ \cdot\ \left({R_o \over {r}}\right)^2$
$\sigma_\theta = P_o + (P_o - \sigma_{re})\ \cdot\ \left({R_o \over {r}}\right)^2$

The displacements in the elastic zone are

$u_r = \left[P_o - \left({2 P_o - q \over {K_p + 1}} \right) \right]\ \left({R_o \over {2G}} \right)\ \left( {R_o \over {r}}\right)$

and in the plastic region are:

$u_r = {r \over {2G}} \chi$
$\chi = (2 \nu - 1)\ \left(P_o + {q \over {K_p -1}}\right)$
$+ {(1 - \nu) (K^2_p - 1) \over {K_p + K_{ps}}} \left(P_i + {q \over {K_p - 1}}\right)\ \left({R_o \over {a}}\right)^{(K_p - 1)}\ \left({R_o \over {r}}\right)^ {(K_{ps} + 1)}$
$+ \left[(1 - \nu)\ {{K_p K_{ps} + 1} \over {K_p + K_{ps}}} - \nu\right]\ \left(P_i + {q \over {K_p - 1}}\right)\ \left({r \over {a}}\right)^{(K_p - 1)}$

where:

$$K_{ps}$$ = $$(1 + \sin \psi)/(1 − \sin \psi)$$ ;

$$\psi$$ = dilation angle; and

$$\nu$$ = Poisson’s ratio.

3DEC Model

The geometry of the 3DEC model used for this analysis is shown in Figure 2. The model is constrained in the direction of the $$x$$-axis, so that it can only move in the $$yz$$-plane (plane-strain condition). The model consists of 900 joined blocks in a radial pattern. The boundary was selected at 10 m (i.e., five hole-diameters) from the hole center. The problem was solved using three different types of zoning. In the first case, a mixed-discretization (hex) zoning was used to discretize 900 blocks into 9000 zones. (An accurate solution of elasto-plastic problems requires finer zoning than is needed for the simulation of linearly elastic problems.) Figure 3 shows the geometry of blocks and zones for this case. (There is one hex zone — i.e., two sets of 5 overlapping tetrahedra per block.) In the second case, 900 blocks were discretized into 50,900 “regular” tetrahedral zones. The geometry of the zones on one side of the model is shown in Figure 4. In all cases, the blocks were joined (i.e., the discontinuities between blocks have no effect on deformation of the model). In the third case, 900 blocks were discretized into 1096 higher-order zones. The geometry of the zones on one side of the model is shown in Figure 5.

Results and Discussion

Because the problem is axisymmetric, all results are presented as a function of the distance from the center of the opening. Figures 6 and 7 compare 3DEC results using hex zoning with the analytical solution. (The radial displacements are compared in Figure 6; the radial and tangential stresses are compared in Figure 7.) The plots show excellent agreement between the two results. Note the 3DEC stresses plotted in Figure 7 are not tetrahedral zone stresses, but stresses obtained by averaging the tetrahedral stresses on the level of the hex zones. (Even in the case of hex zoning, stresses in 3DEC are stored in tetrahedral zones. Quad-zone stresses (i.e., average stresses on the level of a hex zone) are calculated during the calculation of velocities and stress increments only.) A FISH function was used to calculate the average stresses (see “mchole.fis”) .

Results obtained with 3DEC using tetrahedral zoning are compared with the analytical solution in Figures 8 and 9. Agreement between displacements in Figure 8 is excellent. The numerical results deviate from the analytical solution in the vicinity of the far-field boundaries only. (The error is a consequence of the finite size of the numerical model and approximate boundary conditions.) The results for stresses in Figure 9 show large numerical dispersion in the plastic region. Such a result is expected behavior for tetrahedral zones (or triangular zones in 2D) with linear approximation of displacements (or constant approximation of strains and stresses). Mixed discretization, the methodology which is used to improve performance of tetrahedral (or triangular) zones in elastoplastic calculations, is implemented with hex zoning in 3DEC. Results obtained with 3DEC using the high-order tetrahedral zones are compared with the analytical solution in Figures 10 and 11. Agreements between displacements in Figure 10 is excellent. The results for the stresses in Figure 11 are almost as good as the results obtained using the mixed discretization. These results were obtained using only a fraction of the elements needed for the other methods. The different zoning approaches have the following disadvantages:

1. Quad zoning cannot be used for discretization of arbitrarily shaped blocks.

2. Stress fields calculated with tetrahedral zones for elasto-plastic problems are not smooth.

Reference

Salençon, J. “Contraction Quas-Statique D’une Cauite a Symetric Spherique Ou Cylindrique Dans Un Milieu Elastoplastique,” Annales Des Ports Et Chaussees, Vol. 4, pp. 231-236 (1969).

Data Files

mchole-hex.dat

;==========================================================================
; verification test -- 3dec modeling cylindrical hole in an infinite,
;                      isotropic, homogeneous, Mohr-Coulomb medium
; Mohr-Coulomb blocks
;
;==========================================================================
model new

model large-strain off

block create tunnel dip 0 dip-direction 0 radius 1 length 0 0.2 ...
blocks-axial 1 blocks-radial 30 blocks-tangential 15 ...
block delete range cylinder end-1 0 0 0 end-2 0 1 0 radius 0 1
block delete range plane dip 0 dip-direction 0 below
block delete range plane dip 90 dip-direction 90 below

block join on
;
; --- zoning of blocks ---
;
block zone generate hexahedra divisions 1 1 1
;
; --- material properties ---
;
block zone cmodel assign mohr-coulomb
block zone property bulk 3.9e9 shear 2.8e9 density 2600
block zone property cohesion 3.45e6 tension 1e7 friction 30 dilation 30
;
; --- boundary conditions ---
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-x 10
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-z 10
block gridpoint apply velocity-x 0 range position-x 0
block gridpoint apply velocity-y 0 range position-y 0
block gridpoint apply velocity-y 0 range position-y 0.2
block gridpoint apply velocity-z 0 range position-z 0
;
; --- initial conditions ---
;
block insitu stress -30e6 -30e6 -30e6 0 0 0
;
model history mechanical unbalanced-maximum
;
; if there are no joints, local damping should be used
block mechanical damping local 0.8
model solve

program call 'mchole.fis'
model save 'mchole-hex'
;
program return


mchole-tet.dat

;==========================================================================
; verification test -- 3dec modeling cylindrical hole in an infinite,
;                      isotropic, homogeneous, Mohr-Coulomb medium
; tetrahedral zoning
; Mohr-Coulomb blocks
;
;==========================================================================
model new

model large-strain off

model random 10000
block create tunnel dip 0 dip-direction 0 radius 1 length 0 0.2 ...
blocks-axial 1 blocks-radial 10 blocks-tangential 5 ...
block delete range cylinder end-1 0 0 0 end-2 0 1 0 radius 0 1
block delete range plane dip 0 dip-direction 0 below
block delete range plane dip 90 dip-direction 90 below

block join on
;
; --- zoning of blocks ---
;
block zone generate edgelength 0.05 range position-x 0 1.7 position-z 0 1.7
block zone generate edgelength 0.08 range position-x 0 2.9 position-z 0 2.9
block zone generate edgelength 0.15 range position-x 0 5.3 position-z 0 5.3
block zone generate edgelength 0.30
;
; --- material properties ---
;
block zone cmodel assign mohr-coulomb
block zone property bulk 3.9e9 shear 2.8e9 density 2600
block zone property cohesion 3.45e6 tension 1e7 friction 30 dilation 30
;
; --- boundary conditions ---
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-x 10
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-z 10
block gridpoint apply velocity-x 0 range position-x 0
block gridpoint apply velocity-y 0 range position-y 0
block gridpoint apply velocity-y 0 range position-y 0.2
block gridpoint apply velocity-z 0 range position-z 0
;
; --- initial conditions ---
;
block insitu stress -30e6 -30e6 -30e6 0 0 0
;
model history mechanical unbalanced-maximum

; if there are no joints, local damping should be used
block mechanical damping local 0.8
model solve

program call 'mchole.fis'
model save 'mchole-tet'
;
program return


mchole-ho.dat

;==========================================================================
; verification test -- 3dec modeling cylindrical hole in an infinite,
;                      isotropic, homogeneous, Mohr-Coulomb medium
; High order tetrahedral zoning
; Mohr-Coulomb blocks
;
;==========================================================================
model new

model large-strain off

model configure highorder

block create tunnel dip 0 dip-direction 0 radius 1 length 0 0.2 ...
blocks-axial 1 blocks-radial 10 blocks-tangential 5 ...
block delete range cylinder end-1 0 0 0 end-2 0 1 0 radius 0 1
block delete range plane dip 0 dip-direction 0 below
block delete range plane dip 90 dip-direction 90 below

block join on
;
; --- zoning of blocks ---
;
block zone generate edgelength 0.3 range position-x 0.0 1.7 position-z 0.0 1.7
block zone generate edgelength 0.6 range position-x 0.0 2.9 position-z 0.0 2.9
block zone generate edgelength 0.8 range position-x 0.0 5.3 position-z 0.0 5.3
block zone generate edgelength 1
block zone generate high-order-tetra
;
; --- material properties ---
;
block zone cmodel assign mohr-coulomb
block zone property bulk 3.9e9 shear 2.8e9 density 2600
block zone property cohesion 3.45e6 tension 1e7 friction 30 dilation 30
;
; --- boundary conditions ---
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-x 10
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-z 10
block gridpoint apply velocity-x 0 range position-x 0
block gridpoint apply velocity-y 0 range position-y 0
block gridpoint apply velocity-y 0 range position-y 0.2
block gridpoint apply velocity-z 0 range position-z 0
;
; --- initial conditions ---
;
block insitu stress -30e6 -30e6 -30e6 0 0 0
;
model history mechanical unbalanced-maximum
;
; if there are no joints, local damping should be used
block mechanical damping local 0.8
model solve
program call 'mchole.fis'

model save 'mchole-ho'
;
program return


mchole-nmd.dat

;==========================================================================
; verification test -- 3dec modeling cylindrical hole in an infinite,
;                      isotropic, homogeneous, Mohr-Coulomb medium
; Nodal Mixed Discretization zoning
; Mohr-Coulomb blocks
;
;==========================================================================
model new

model large-strain off

block create tunnel dip 0 dip-direction 0 radius 1 length 0 0.2 ...
blocks-axial 1 blocks-radial 10 blocks-tangential 5 ...
block delete range cylinder end-1 0 0 0 end-2 0 1 0 radius 0 1
block delete range plane dip 0 dip-direction 0 below
block delete range plane dip 90 dip-direction 90 below

block join on
;
; --- zoning of blocks ---
;
block zone generate edgelength 0.1 range position-x 0 1.7 position-z 0 1.7
block zone generate edgelength 0.16 range position-x 0 2.9 position-z 0 2.9
block zone generate edgelength 0.3 range position-x 0 5.3 position-z 0 5.3
block zone generate edgelength 0.6
;
; --- material properties ---
;
block zone cmodel assign mohr-coulomb
block zone property bulk 3.9e9 shear 2.8e9 density 2600
block zone property cohesion 3.45e6 tension 1e7 friction 30 dilation 30

block zone nodal-mixed-discretization on
;
; --- boundary conditions ---
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-x 10
block face apply stress -30e6 -30e6 -30e6 0 0 0 range position-z 10
block gridpoint apply velocity-x 0 range position-x 0
block gridpoint apply velocity-y 0 range position-y 0
block gridpoint apply velocity-y 0 range position-y 0.2
block gridpoint apply velocity-z 0 range position-z 0
;
; --- initial conditions ---
;
block insitu stress -30e6 -30e6 -30e6 0 0 0
;
model history mechanical unbalanced-maximum

; if there are no joints, local damping should be used
block mechanical damping local 0.8
model solve

program call 'mchole.fis'

model save 'mchole-nmd'

program return


mchole.fis

;---------------------------------------------------------------------
;    exact and numerical solution processing for salencon problem
;
;    nastr, stores in
;    Table 1: analytical values rad/a -sigr/p0
;    Table 2: analytical values rad/a -sigt/p0
;    Table 3: numerical values rad/a -sigr/p0
;             at zone centroid closest to x axis
;    Table 4: numerical values rad/a -sigt/p0
;             at zone centroid closest to x axis
;
;    and calculates in whole grid
;    errsr  : average %error in -sigr/p0
;    errst  : average %error in -sigt/p0
;
;    Table 5: analytical values of rad/a -ur/a
;    Table 6: numerical  values of rad/a -ur/a at grid points on x axis
;
;    and calculates in whole grid
;    errd   : average %error in modulus of displacement
;---------------------------------------------------------------------
fish automatic-create off
;
;--- access property (cohesion, friction) and grid information ---
fish define const
global po = 30.0e6 ; far-field stress
global kp
global kpm1
global cqk
global sre
global ro
global kpsp1
global xi1
global xi2
global xi3
global coe

; material properties
local cco = 3.45e6
local bm    = 3.3e9
local sm    = 2.8e9
local nu    = (3.*bm-2.*sm)/(6.*bm+2.*sm)

kp    = (1.+math.sin(cfi))/(1.-math.sin(cfi))
kpm1  = kp - 1
local cq    = 2.*cco*math.cos(cfi)/((1.0-math.sin(cfi))*po)
cqk   = cq/(kp-1.)
sre   = (2.-cq)/(kp+1.)
ro    = (2.*(1.+cqk)/((kp+1.)*cqk))^(1./kpm1)
local kps   = (1.+math.sin(cdi))/(1.-math.sin(cdi))
kpsp1 = kps+1
xi1   = (2.*nu-1.)*(1.+cqk)
xi2   = ((1.-nu)*(kp^2-1.)/(kp+kps))*cqk*ro^kpm1
xi3   = ((1.-nu)*(kp*kps+1)/(kp+kps)-nu)*cqk
coe   =  po/(2.*sm)
end
; =================================================
; Stresses
; =================================================
fish define nastr
local tab1  = table.get(1)
local tab2  = table.get(2)
local tab3  = table.get(3)
local tab4  = table.get(4)
table.label(tab1) = 'analytic-sigr'
table.label(tab2) = 'analytic-sigt'
table.label(tab3) = '3DEC-sigr'
table.label(tab4) = '3DEC-sigt'
global errsr
global errst
errsr = 0.0
errst = 0.0

loop foreach local iz_ block.zone.list
local mark = 0
local igp
loop igp (1,4)
local gpind = block.zone.gp(iz_,igp)
if block.gp.pos.z(gpind) < 0.001 then
mark = mark + 1
end_if
end_loop
local rad = math.sqrt(block.zone.pos.x(iz_)^2 + ...
block.zone.pos.z(iz_)^2)
local theta_ = math.atan2(block.zone.pos.z(iz_),block.zone.pos.x(iz_))
local sigre = -cqk*(1.-aux1)
local sigte = -cqk*(1.-kp*aux1)
else
sigre = 1. - aux1
sigte = 1. + aux1
end_if
local stemp1 = 0.5*(block.zone.stress.xx(iz_) + ...
block.zone.stress.zz(iz_))
local stemp2 = 0.5*(block.zone.stress.xx(iz_) - ...
block.zone.stress.zz(iz_))
local stemp3 = block.zone.stress.xz(iz_)*math.sin(2.0*theta_)
local sigr = -(stemp1 + stemp2*math.cos(2.0*theta_)+stemp3)/po
local sigt = -(stemp1 - stemp2*math.cos(2.0*theta_)-stemp3)/po
if mark > 2
end_if
local err = math.abs(sigr - sigre)
errsr = errsr + err
err = math.abs(sigt - sigte)
errst = errst + err
end_loop
errsr = errsr * 100. / block.zone.num
errst = errst * 100. / block.zone.num
end
; =================================================
; Displacements
; =================================================
local tab5  = table.get(5)
local tab6  = table.get(6)
table.label(tab5) = 'analytic-ur'
table.label(tab6) = '3DEC-ur'
global errd
errd = 0.0

loop foreach local igp_ block.gp.list
local mark = 0
if block.gp.pos.z(igp_) < 0.001 then
if block.gp.pos.y(igp_) < 0.001 then
mark = 1
end_if
end_if
else
end_if
local dis = math.sqrt(block.gp.disp.x(igp_)^2 + ...
block.gp.disp.z(igp_)^2)
local err = math.abs(dis - dise)
errd = errd + err
if mark = 1 then
end_if
end_loop
errd = errd * 100. / (dism * block.zone.num)
end
;
[const]
[nastr]