---
jupyter:
  description: 'Stress in a thick-walled pipe with planex2d and NGSolve: Lamé''s radial and hoop stress in a quarter of
    the pipe, the bore stress that never falls below the pressure however thick the wall, and a shrink-fitted ring,
    checked against the compound-cylinder formulas, that lowers the peak stress by a fifth.'
---

# Stress in a thick-walled pipe

Why does a thicker wall help a pressure pipe so little? Lamé's solution for a thick-walled cylinder gives the hoop
stress across a wall from the bore radius $a$ to the outer radius $b$ under the pressure $p$ as

$$\sigma_\theta(r) = \frac{p a^2}{b^2 - a^2} \left(1 + \frac{b^2}{r^2}\right), \qquad
\sigma_\theta(a) = p \, \frac{b^2 + a^2}{b^2 - a^2}.$$

The stress is highest at the bore, and however thick the wall, it never falls below the pressure itself, because the
outer layers carry ever less. The mean over the wall is $p a / (b - a)$ for any thickness, as the balance of forces on
half the pipe demands, so the thin-wall formula is right on average and wrong at the bore. Gun barrels and high-pressure
vessels get around this with layers shrunk onto each other, which put the bore in compression before the pressure comes.

The example is a steel pipe of bore radius 50 and outer radius 100, in millimeters, under 100 MPa. NGSolve solves a
quarter of it, and its stresses are checked against Lamé's solution. The bore stress is then computed for walls of
growing thickness, and finally the same wall is made of two layers, the outer one shrunk onto the inner, and checked
against the formulas for a compound cylinder.

This page is part of a collection of [physics examples](physics.md) using planex2d.

![The hoop stress in a pipe made of a tube with a ring shrunk onto it, under pressure](../_static/gallery/thick-walled-pipe.png){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
```

```python nbsphinx="hidden"
from docs_figures import MAP, Chart, Row
```

## The pipe

A quarter of the pipe is enough, because its straight sides lie on lines of symmetry, where the displacement across them
is zero. The wall can be made of layers, named `pipe` for the first and `ring` for a second one:

```python
E, nu = 210e3, 0.3
a, b, p = 50, 100, 100


def quarter(r0, r1, name):
    """a quarter of the ring between the radii r0 and r1"""
    sk = px.QuickSketch((r0, 0)).Line(r1 - r0)
    sk.Arc(center=(0, 0), sweep=90).Turn(90).Line(r1 - r0)
    face = sk.Arc(center=(0, 0), sweep=-90).Face()
    face.name = name
    return face


def pipe(*radii):
    """a quarter of the pipe with its layers between the radii given, its bore and symmetry edges named"""
    shape = px.Glue([quarter(r0, r1, name) for r0, r1, name in zip(radii, radii[1:], ("pipe", "ring"))])
    for edge in shape.edges:
        (x0, y0), (x1, y1) = edge.bounding_box
        if y1 < 1e-9:
            edge.name = "bottom"
        elif x1 < 1e-9:
            edge.name = "left"
        elif max(x1, y1) < radii[0] + 1e-9:
            edge.name = "bore"
    return shape


single = pipe(a, b)
# px.Plot(single)
```

The wall is in plane stress, with the pressure on the bore. A ring that is shrunk on is made too small by the relative
amount `shrink`, and its stress comes only from the strain it has beyond that one. The function returns the radial and
the hoop stress:

```python
mu, lam = E / (2 * (1 + nu)), E * nu / (1 - nu**2)


def stresses(shape, pressure, shrink=0.0, order=3):
    """the mesh and the radial and hoop stress, with the pressure on the bore and the ring made too small by shrink"""
    mesh = ngs.Mesh(px.OCCGeometry(shape).GenerateMesh(maxh=5))
    mesh.Curve(order)
    V = ngs.VectorH1(mesh, order=order, dirichletx="left", dirichlety="bottom")
    u, v = V.TnT()

    def stress(strain):
        return 2 * mu * strain + lam * ngs.Trace(strain) * ngs.Id(2)

    def strain(w):
        return ngs.Sym(ngs.Grad(w))

    made = (mesh.MaterialCF({"ring": -shrink}, default=0) if shrink else ngs.CF(0)) * ngs.Id(2)
    A = ngs.BilinearForm(ngs.InnerProduct(stress(strain(u)), strain(v)) * ngs.dx).Assemble()
    f = ngs.LinearForm(V)
    if pressure:
        f += -pressure * ngs.specialcf.normal(2) * v * ngs.ds("bore")
    if shrink:
        f += ngs.InnerProduct(stress(made), strain(v)) * ngs.dx("ring")
    f.Assemble()
    w = ngs.GridFunction(V)
    w.vec.data = A.mat.Inverse(V.FreeDofs(), inverse="sparsecholesky") * f.vec
    s = stress(strain(w) - made)
    c, d = ngs.x / ngs.sqrt(ngs.x**2 + ngs.y**2), ngs.y / ngs.sqrt(ngs.x**2 + ngs.y**2)
    radial = s[0, 0] * c * c + s[1, 1] * d * d + 2 * s[0, 1] * c * d
    hoop = s[0, 0] * d * d + s[1, 1] * c * c - 2 * s[0, 1] * c * d
    return mesh, radial, hoop


def at(mesh, cf, r, angle=30):
    """a stress at the radius r, on the ray at the angle given, just inside the wall at its bore and outside"""
    r = min(max(r, a + 1e-6), b - 1e-6)
    t = math.radians(angle)
    return cf(mesh(r * math.cos(t), r * math.sin(t)))
```

The computed stresses follow Lamé's solution across the wall, and the hoop stress integrated over the wall gives
$p a$, the force per unit length that half the pipe has to carry:

```python
single_mesh, single_radial, single_hoop = stresses(single, p)
k = p * a**2 / (b**2 - a**2)
for r in (a, 60, 75, 90, b):
    print("r", r, "| radial", at(single_mesh, single_radial, r), "Lamé", k * (1 - b**2 / r**2),
          "| hoop", at(single_mesh, single_hoop, r), "Lamé", k * (1 + b**2 / r**2))
rs = np.linspace(a, b, 2001)
print("hoop stress over the wall", np.trapezoid([at(single_mesh, single_hoop, r) for r in rs], rs), "| p a", p * a)
# px.PlotField(single_hoop, single_mesh, single)
```

## The wall thickness

With the outer radius from 1.1 to 5 times the bore radius, the stress at the bore falls towards the pressure, 100 MPa,
and no further. The thin-wall formula $p a / (b - a)$ falls to zero (Figure 1):

```python
ratios = [1.1, 1.25, 1.5, 2, 3, 4, 5]
bore = []
for ratio in ratios:
    m, _, h = stresses(pipe(a, ratio * a), p)
    bore.append(at(m, h, a))
    print("b / a", ratio, "| hoop stress at the bore", bore[-1], "| Lamé", p * (ratio**2 + 1) / (ratio**2 - 1),
          "| thin wall", p / (ratio - 1))
# plt.plot(ratios, bore); plt.show()
```

```python nbsphinx="hide-input"
fine = np.linspace(1.1, 5, 200)
Chart("Outer radius over bore radius", "Hoop stress at the bore in MPa",
      lines={"Lamé": (fine, p * (fine**2 + 1) / (fine**2 - 1)), "Thin wall": (fine, p / (fine - 1)),
             "The pressure": (fine, np.full_like(fine, p))},
      dots={"NGSolve": (ratios, bore)}, title="The bore stress against the wall thickness")
```

## A shrink-fitted ring

The same wall is now made of two layers, a tube from 50 to 75 and a ring from 75 to 100. The ring is made a little too
small and shrunk on, so it presses on the tube with the pressure $p_c$, which for an interference $\delta$ at the radius
$c$ between the layers is

$$p_c = \frac{E \delta}{2 c^3} \, \frac{(b^2 - c^2)(c^2 - a^2)}{b^2 - a^2}.$$

It puts the tube's bore in compression, by $2 p_c c^2 / (c^2 - a^2)$, and stretches the ring's bore by
$p_c (b^2 + c^2) / (b^2 - c^2)$. Under the pressure, both add to Lamé's stresses of the whole wall, and the
interference is chosen from these formulas so that the two bores carry the same hoop stress:

```python
c = 75
tube_bore, ring_bore = k * (1 + b**2 / a**2), k * (1 + b**2 / c**2)
pc = (tube_bore - ring_bore) / (2 * c**2 / (c**2 - a**2) + (b**2 + c**2) / (b**2 - c**2))
delta = 2 * c**3 * (b**2 - a**2) * pc / (E * (b**2 - c**2) * (c**2 - a**2))
print("interference", delta, "| contact pressure", pc,
      "| hoop stress at both bores", tube_bore - 2 * pc * c**2 / (c**2 - a**2))
```

Shrunk on alone, the computed ring presses on the tube with the contact pressure of the formula. Under the pressure,
the peak hoop stress, 167 MPa in the single wall, drops by about a fifth and sits at both bores (Figure 2):

```python
layered = pipe(a, c, b)
mesh, radial, hoop = stresses(layered, 0, shrink=delta / c)
print("shrunk on: contact pressure", -at(mesh, radial, c), "| formula", pc,
      "| tube's bore", at(mesh, hoop, a), "| formula", -2 * pc * c**2 / (c**2 - a**2))
mesh2, radial2, hoop2 = stresses(layered, p, shrink=delta / c)
print("under pressure: tube's bore", at(mesh2, hoop2, a), "| ring's bore", at(mesh2, hoop2, c + 1e-3),
      "| single wall's bore", tube_bore)
# px.PlotField(hoop2, mesh2, layered)
```

```python nbsphinx="hide-input"
def sampled(mesh, cf):
    """the hoop stress for Row, at points of the face, read in the quarter that was solved"""
    def field(face, plot=True, **_):
        (x0, y0), (x1, y1) = face.bounding_box
        gx, gy = np.meshgrid(np.linspace(x0, x1, 400), np.linspace(y0, y1, 400))
        xy = np.c_[gx.ravel(), gy.ravel()]
        xy = xy[face.Distance(xy, signed=True) < -0.05]
        values = np.asarray(cf(mesh(np.abs(xy[:, 0]), np.abs(xy[:, 1])))).ravel()
        return {"points": [tuple(q) for q in xy], "values": values, "label": "hoop stress in MPa"}
    return field


Row((single, "One wall", {"field": sampled(single_mesh, single_hoop), "clim": (60, 170), "cmap": MAP.reversed()}),
    (layered, "Tube and shrunk-on ring", {"field": sampled(mesh2, hoop2), "clim": (60, 170), "cmap": MAP.reversed()}),
    caption="The hoop stress under 100 MPa, in one wall and in two layers")
```

Across the wall, the single wall's hoop stress falls from the bore outwards. The layered wall's starts lower, jumps up
at the ring's bore and falls again (Figure 3):

```python
radii = np.linspace(a, b, 41)
across = {"One wall": [at(single_mesh, single_hoop, r) for r in radii],
          "Tube and ring": [at(mesh2, hoop2, r + (1e-3 if r == c else 0)) for r in radii]}
# plt.plot(radii, across["Tube and ring"]); plt.show()
```

```python nbsphinx="hide-input"
inner, outer = np.linspace(a, c, 100), np.linspace(c, b, 100)
layered_lame = np.r_[k * (1 + b**2 / inner**2) - pc * c**2 / (c**2 - a**2) * (1 + a**2 / inner**2),
                     k * (1 + b**2 / outer**2) + pc * c**2 / (b**2 - c**2) * (1 + b**2 / outer**2)]
Chart("Radius in mm", "Hoop stress in MPa",
      lines={"One wall, Lamé": (np.r_[inner, outer], k * (1 + b**2 / np.r_[inner, outer] ** 2)),
             "Tube and ring, Lamé": (np.r_[inner, outer], layered_lame)},
      dots={"NGSolve": (np.r_[radii, radii], np.r_[across["One wall"], across["Tube and ring"]])},
      title="The hoop stress across the wall")
```

See [Stress at a shoulder fillet](shoulder-fillet.md) for a stress concentration, and
[Cantilever deflection](cantilever-deflection.md) for a plane solid against beam theory.

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from docs_figures import INK
from gallery import save_figure

full = px.Glue([px.Circle((0, 0), c).Face() - px.Circle((0, 0), a).Face(),
                px.Circle((0, 0), b).Face() - px.Circle((0, 0), c).Face()])
fig, ax = plt.subplots(figsize=(6, 4))
px.PlotField(sampled(mesh2, hoop2), full, ax=ax, colormap=MAP.reversed(), min=60, max=170)
full.Plot(ax=ax, face_col=None, line_col=INK, arc_col=INK, legend=False)
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(-157, 157)
ax.set_ylim(-105, 105)
save_figure(fig, "thick-walled-pipe", bare=True, bars={1: "hoop stress in MPa"})
plt.close(fig)
```