Hoek-Brown tunnel: the closed-form ground reaction curve

There is a kind of confidence you can only get from a problem that has an exact answer. Most geotechnical finite element work does not: you mesh a real slope, you run it, and the number that comes out is the only number you have. Nothing tells you whether it is right. So the discipline is to step away from the real problems occasionally and solve one that a mathematician solved first, where disagreement is not a matter of opinion.

This is one of those. A circular tunnel is driven in a rock mass under uniform in-situ stress, and the rock is described by the Hoek-Brown failure criterion. The problem has a closed-form elasto-plastic solution. Below it is derived, the same problem is built in KATAI 2D through its Python API, and the two are put side by side. KATAI 2D is software I write and maintain, so read this as a documented verification record rather than an independent review.

The short version: the finite element radial stress lands on the closed form to +0.04% at the elastic-plastic boundary and +0.18% in the far field, and the answer moves by 0.04% when the mesh is refined by half. The longer version is more interesting, because building this comparison found three real defects in the program — one of which had been invisible to three separate verification tests. They are described below and all three are fixed, which is the entire reason for solving a problem whose answer is already known.

Why is rock not soil?Permalink to “Why is rock not soil?”#

Almost every constitutive model a geotechnical engineer uses day to day describes strength with a straight line: Mohr-Coulomb, with a cohesion $c'$ and a friction angle $\varphi'$. For soil that is a good description. For rock it is not, and the reason is worth stating precisely rather than as a slogan.

Rock strength is strongly non-linear in the confining stress. Squeeze a rock sample harder from the sides and it does not get stronger in proportion — it gets stronger more slowly, and the strength envelope curves over. The Hoek-Brown criterion (Hoek, Carranza-Torres and Corkum, 2002) describes that curve directly:

$\sigma_1 = \sigma_3 + \sigma_{ci} \left( m_b \frac{\sigma_3}{\sigma_{ci}} + s \right)^{a}$

with compression positive and $\sigma_1 \geq \sigma_3$, which is the rock-mechanics convention and the one used throughout this article. A finite element program normally works in the opposite one, compression negative, where the same criterion reads $\sigma_1 = \sigma_3 - |\sigma_{ci}|\,(m_b(-\sigma_3)/|\sigma_{ci}| + s)^a$ with the roles of the two principal stresses exchanged. The physics is identical; only the bookkeeping differs, and mixing the two in one text is a good way to lose an afternoon.

Here $\sigma_{ci}$ is the uni-axial compressive strength of the intact rock, and $m_b$, $s$, $a$ are rock-mass constants that follow from the Geological Strength Index $GSI$ and a disturbance factor $D$:

Fitting a straight line to that curve is not an approximation that is merely a bit conservative. It is right at the confining stress you fitted it at, and wrong on both sides of it — unconservative in one direction and over-conservative in the other, and which direction depends on where you fitted. That is not a comfortable position for a design.

Note also what these parameters are. $\sigma_{ci}$ comes from a laboratory test on a core. $m_i$ is a property of the rock type. $GSI$ is a field judgement about how blocky the mass is and how weathered its joint surfaces are. $D$ is a judgement about how much the excavation damaged it. None of them is a fitted $c'$–$\varphi'$ pair, and that is the point: the model asks the questions a geologist can actually answer.

What is the problem, and where do its numbers come from?Permalink to “What is the problem, and where do its numbers come from?”#

A tunnel of radius $R$ is excavated in a rock mass under a hydrostatic in-situ stress $\sigma_0$. Before excavation the rock is uniformly compressed. As the tunnel is driven and the support pressure $p_i$ on the wall falls from $\sigma_0$ towards zero, the rock relaxes inward. Below a critical support pressure a plastic annulus opens around the opening and grows as $p_i$ falls further.

Plotting the wall displacement against $p_i$ gives the ground reaction curve, the central object of the convergence-confinement method of tunnel design. The elasto-plastic solution for a Hoek-Brown medium is due to Carranza-Torres and Fairhurst (1999, 2000), with the generalised form ($a \neq \tfrac12$) in Carranza-Torres (2004).

The rock mass and stress state used here are case D₁ of Table A1.1 in Hoek, Carranza-Torres, Diederichs and Corkum (2008), Integration of geotechnical and structural design in tunnelling, 56th Annual Geotechnical Engineering Conference, University of Minnesota. That case is a graphitic phyllite at about 1100 m depth, from the Yacambú-Quibor tunnel in Venezuela — a genuinely difficult piece of ground.

The same table publishes the derived rock-mass constants, which gives a first check before anything is solved at all. Computing them from $GSI$ and $D$:

Agreement to the printed digits. The arithmetic is the same arithmetic.

What is the closed-form solution?Permalink to “What is the closed-form solution?”#

I derive it here rather than quote it, for a reason that matters more than it might seem: a formula copied from memory into a verification is not a verification, it is a comparison of a program against a recollection.

Compression positive. Radial equilibrium of an axisymmetric ring, with no body force, is

$\frac{d\sigma_r}{dr} = \frac{\sigma_\theta - \sigma_r}{r}$

Inside the plastic annulus the stress pair sits on the Hoek-Brown envelope, so

$\sigma_\theta = \sigma_r + \sigma_{ci}\left(\frac{m_b \sigma_r}{\sigma_{ci}} + s\right)^{a}$

Substituting the second into the first leaves an equation in $\sigma_r$ and $r$ alone, and the variables separate:

$\frac{d\sigma_r}{\sigma_{ci}\left(\dfrac{m_b\sigma_r}{\sigma_{ci}} + s\right)^{a}} = \frac{dr}{r}$

With the substitution $u = m_b\sigma_r/\sigma_{ci} + s$ this integrates in closed form for any $a \neq 1$. Applying the boundary condition $\sigma_r = p_i$ at $r = R$:

$\ln\frac{r}{R} = \frac{1}{m_b(1-a)}\left[\left(\frac{m_b\sigma_r}{\sigma_{ci}} + s\right)^{1-a} - \left(\frac{m_b p_i}{\sigma_{ci}} + s\right)^{1-a}\right]$

That is the whole plastic branch. What remains is to find where it ends.

Outside the plastic annulus the rock is elastic, and the elastic field around a hole in an infinite medium under hydrostatic $\sigma_0$ has the property that the two in-plane normal stresses always sum to $2\sigma_0$. At the elastic-plastic boundary the stress must satisfy that and lie on the envelope, so the radial stress there, $\sigma_r^{ep}$, solves

$2\left(\sigma_0 - \sigma_r^{ep}\right) = \sigma_{ci}\left(\frac{m_b\sigma_r^{ep}}{\sigma_{ci}} + s\right)^{a}$

Put $\sigma_r^{ep}$ into the integral above and you have the plastic radius $R_{pl}$. Note what this second equation also gives: at $p_i = \sigma_r^{ep}$ the plastic radius equals $R$, so $\sigma_r^{ep}$ is the critical support pressure below which yielding starts at all.

How do we know the derivation is right?Permalink to “How do we know the derivation is right?”#

A derivation checked only against itself is worth very little. Set $a = \tfrac12$ — the original Hoek-Brown criterion — and introduce the scaled stresses Carranza-Torres and Fairhurst use,

$S = \frac{\sigma}{m_b\sigma_{ci}} + \frac{s}{m_b^2}$

Then $\sqrt{S} = \frac{1}{m_b}\sqrt{m_b\sigma/\sigma_{ci} + s}$, and the integral collapses to

$R_{pl} = R\,\exp\!\left[2\left(\sqrt{S^{ep}} - \sqrt{P_i}\right)\right]$

while the elastic-matching equation becomes a quadratic in $\sqrt{S}$ whose root is

$P_i^{cr} = \frac{1}{16}\left[1 - \sqrt{1 + 16 S_0}\right]^2$

These are exactly the published scaled forms. Two independent routes — a fresh integration here and the scaled solution in the literature — arrive at the same expressions, so neither is a transcription of the other.

For this case the numbers are a critical support pressure of 14.139 MPa, which is 50.5% of the in-situ stress — this rock starts yielding while the tunnel is still heavily supported — and a plastic radius when fully unsupported of 2.552 R.

How do you model a circular tunnel as a one-dimensional strip?Permalink to “How do you model a circular tunnel as a one-dimensional strip?”#

The obvious way to model a circular tunnel is to draw a circle, mesh a quarter of it with symmetry conditions, and apply pressure to the arc. That works, but it means approximating the arc with facets and applying one load per facet, each with its own normal direction.

There is a cleaner route, and it is worth explaining because it looks at first like a different problem entirely. In an axisymmetric analysis the three normal stresses are $(\sigma_r, \sigma_z, \sigma_\theta)$ — radial, axial and hoop. Map the tunnel onto that: let $r$ be the tunnel radius, $\theta$ the hoop direction, and $z$ the tunnel axis. Holding $u_z = 0$ on both horizontal edges is then precisely the tunnel’s plane-strain condition along its own axis. The whole cross-section reduces to a one-dimensional radial strip with straight edges: inner face at $r = R$, outer face at $r = B$, two pressure boundaries instead of forty facet loads.

import katai

prj = katai.Project("Rock tunnel", kind="axisymmetric",
                    mesh_size=0.4, auto_refine=False)

rock = prj.materials.hoek_brown("Graphitic phyllite", E=7.5e6, nu=0.25,
                                sigma_ci=50.0e3, mi=7.0, gsi=48.0, D=0.0,
                                psi=0.0, sigma_psi=0.0, gamma=0.001)

# r = R is the tunnel wall, r = B the far field; u_z = 0 top and bottom
# is plane strain along the tunnel axis.
prj.geometry.rectangle(2.5, 0.0, 100.0, 2.0, material=rock,
                       bottom="vertical", top="vertical",
                       left="free", right="free")

prj.loads.line_load((100.0, 0.0), (100.0, 2.0), qx=-28.0e3,
                    name="In-situ stress")

The support pressure is then let down in stages. Each stage is one load and one phase, so that removing it is a plain deactivation and the phase’s own displacement is that stage’s convergence:

pieces, previous = [], 28.0e3
for target in LEVELS:            # [24500, 21000, ..., 60, 30, 0] kPa
    pieces.append(prj.loads.line_load((2.5, 0.0), (2.5, 2.0),
                                      qx=+(previous - target)))
    previous = target

for piece, target in zip(pieces, LEVELS):
    prj.phases.plastic(f"p_i = {target/1000:.2f} MPa", deactivate=[piece],
                       steps=90, max_iterations=1500)

That is the whole model, closed form included, in about thirty lines.

Why verify the fixture before trusting it?Permalink to “Why verify the fixture before trusting it?”#

Here is a habit worth borrowing. Before putting the rock model into this fixture, run it with a material whose answer is already known exactly, and check that the fixture is right. Otherwise a disagreement later has two possible causes and you cannot tell them apart.

So: same geometry, same boundaries, same axisymmetric mapping, but a linear elastic material, checked against the Lamé thick-cylinder solution. Result: finite element $-1.170768 \times 10^{-2}$ m against Lamé’s $1.171053 \times 10^{-2}$ m — a difference of 0.024%. The fixture is right, and any later disagreement belongs to the rock model rather than to the geometry.

What does the comparison give?Permalink to “What does the comparison give?”#

The ground reaction curvePermalink to “The ground reaction curve”#

Nineteen converged stages, from a nearly fully supported tunnel down to a support pressure of 0.06 MPa, which is 0.2% of the in-situ stress:

The elastic part of that curve is itself a check. While the rock is still elastic each 3.5 MPa stage should close the wall by exactly

$\Delta u = \frac{(1+\nu)\,\Delta p\,R}{E} = \frac{1.25 \times 3500 \times 2.5}{7.5\times10^6} = 1.4583\ \text{mm}$

and the first four stages give 1.46 mm each. Then, at 14 MPa — the critical pressure the closed form predicted — the curve begins to bend away from the straight line. The program starts yielding where the mathematics says it should.

The radial stress profilePermalink to “The radial stress profile”#

This is the sharper test, because it compares a whole field rather than a single number. At $p_i = 0.06$ MPa:

The row in bold is the elastic-plastic boundary. The finite element solution passes through it at +0.04%, having come through the plastic annulus on one side and the elastic field on the other, with no knowledge of where the boundary was supposed to be. It found it.

One honest note: at the tunnel wall itself the recovered nodal stress is 0.054 MPa against the applied 0.060 MPa — a gap of 0.006 MPa, which is 0.02% of the in-situ stress. The program prints that gap as −10.50%, and the percentage is the misleading way to read it: it is a large fraction of a very small number. Nodal stress extrapolation at a free boundary is the weakest point of any displacement-based finite element solution, and it is worth quoting in the unit that shows what it actually is.

Does the answer move with the mesh?Permalink to “Does the answer move with the mesh?”#

A number that changes when you refine the mesh is not an answer, it is a coincidence. Halving the element size from 0.6 m to 0.4 m:

The wall closure is stable to 0.04%. The plastic-annulus stress error shrinks with refinement, which is what a discretisation error does and a model error does not.

How does it compare with the published finite element results?Permalink to “How does it compare with the published finite element results?”#

The same source publishes its own finite element results for this case (Table A1.3): a normalised plastic radius of 2.3 and a maximum wall displacement of 49.5 mm.

The displacements agree to about 2%, which for two independently written programs on an elasto-plastic problem is close.

The plastic radius deserves a comment rather than a shrug. The closed form gives 2.552, and that value is exact for the stated parameters — the plastic radius is a property of the stress solution alone and does not depend on the flow rule, the dilatancy or the mesh. The published 2.3 corresponds to a support pressure of about 0.4 MPa on the closed-form curve. The most likely explanation is how a plastic radius gets read off a finite element mesh: yield indicators are discrete, and the outermost ring of slightly yielding elements is easy to miss. This is a good illustration of why a closed form is worth the trouble. It is the reference against which both programs are measured, rather than one program being measured against the other.

What did the comparison find that the element tests missed?Permalink to “What did the comparison find that the element tests missed?”#

Now the part that is uncomfortable to write, and the reason the exercise was worth doing at all. Building this comparison found three defects in the program, and they are all instances of one pattern: making a model reachable turns every place its name is not written into a bug.

One — the model was never integrated in axisymmetry at all. The constitutive dispatch for axisymmetric analysis is a switch over a closed enumeration. The Hoek-Brown branch had simply never been written, and the switch had no default: case. So for a rock material in axisymmetry the trial stress and the tangent were left exactly as the caller passed them in — the previous iterate, handed back as though it had been integrated. Nothing failed loudly. The equilibrium iteration merely did not converge, and the run reported a stalled Newton search, which is also what a genuinely difficult problem looks like.

What makes this one instructive is that the model had already been verified three times: at the material point against closed forms, and through the whole finite element path from an input file. All three of those tests are plane strain, and all three are element tests. A missing branch in a switch over a closed enum is about the cheapest silent failure there is, and it survived a verification regime that was otherwise fairly thorough.

The fix is one branch. The lasting fix is a gate that walks the model registry and pushes every registered model through both integrators, with the outputs poisoned with NaN first — so a branch that does not exist cannot pass by handing back the caller’s own values. The gate was checked by removing the new branch again and confirming that it goes red.

Two — the tangent was not strong enough for a boundary value problem. The rock model handed the solver the elastic operator rather than a consistent tangent. That is enough to integrate the material correctly, and every material-point test passes with it, because those tests ask what stress comes back. It is not enough to iterate an ill-conditioned boundary value problem to equilibrium: the tunnel stalled as soon as the plastic annulus passed about 4% of the radius, and more load steps did not help — 150 steps per stage reached exactly the same wall as 40. Newton with a wrong-but-symmetric operator converges linearly, and linearly is not always convergent. The fix was a finite-difference consistent tangent, which two other models in the same program already use.

Three — and this one is the sharpest — the retry that uses that tangent never reached the rock. The solver has a hybrid strategy: when an increment fails with the cheap tangent, retry that increment with the stronger one. The flag that arms it was written when one model was the only one with a stronger tangent, and it tested for that model by name. So the new tangent existed and was never called.

The evidence for this was itself worth the price of admission. After adding the tangent, the run stalled at the identical load factor as before — 0.750 — and at 40 steps as at 150. By the solver’s own diagnostic reasoning, a limit load does not move with the load-step count, and a stall does; a number that does not move at all meant the new code was not on the path. It was not.

The flag’s name was the tell. It asked “is Hardening Soil present?” when the question it needed to ask was “is there a material here that can be asked for a stronger tangent?” A predicate that tests for an instance where it means a property will be wrong the first time somebody adds a second instance — and it will be wrong silently.

With all three fixed, the run goes from 0% of the unloading to nineteen converged stages and the agreement in the tables above. All three are now pinned by a regression case that asserts two things deliberately: that the unloading converges, and that where it converges it is right. Either one alone would have passed at some point during the repair.

Reproducing itPermalink to “Reproducing it”#

The full script ships with the program as python/examples/rock_tunnel_ground_reaction.py, closed form included, and it prints every table above. It needs no data files and asserts its own results, so it fails the test suite rather than the first reader who tries it.

The last 0.2% of the unloading — from 0.06 MPa to exactly zero — still does not converge, and the run says so rather than reporting a partial result as an answer. That is the next thing to fix, and it will be measured the same way.

ReferencesPermalink to “References”#

• Hoek, E., Carranza-Torres, C. and Corkum, B. (2002). Hoek-Brown failure criterion — 2002 edition. Proceedings of NARMS-TAC, Toronto, 267–273.
• Carranza-Torres, C. and Fairhurst, C. (1999). The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Sciences 36(6), 777–809.
• Carranza-Torres, C. and Fairhurst, C. (2000). Application of the convergence-confinement method of tunnel design to rock masses that satisfy the Hoek-Brown failure criterion. Tunnelling and Underground Space Technology 15(2), 187–213.
• Carranza-Torres, C. (2004). Elasto-plastic solution of tunnel problems using the generalized form of the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Sciences 41(S1), 629–639.
• Hoek, E., Carranza-Torres, C., Diederichs, M. and Corkum, B. (2008). Integration of geotechnical and structural design in tunnelling. 56th Annual Geotechnical Engineering Conference, University of Minnesota.