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 using planex2d.

import math

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

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:

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:

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:

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)

r 50 | radial -99.99247320096124 Lamé -100.0 | hoop 166.66101918033394 Lamé 166.66666666666669
r 60 | radial -59.256076509991836 Lamé -59.25925925925926 | hoop 125.90918864357161 Lamé 125.92592592592594
r 75 | radial -25.927000822138098 Lamé -25.925925925925924 | hoop 92.59472965111968 Lamé 92.5925925925926
r 90 | radial -7.818762292323026 Lamé -7.818930041152262 | hoop 74.48648910942858 Lamé 74.48559670781894
r 100 | radial 0.0008809149583157705 Lamé 0.0 | hoop 66.66682818958583 Lamé 66.66666666666667
hoop stress over the wall 5000.000415239379 | p a 5000

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):

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()

b / a 1.1 | hoop stress at the bore 1052.3746115619701 | Lamé 1052.3809523809514 | thin wall 999.9999999999991
b / a 1.25 | hoop stress at the bore 455.5555502238628 | Lamé 455.55555555555554 | thin wall 400.0
b / a 1.5 | hoop stress at the bore 259.9939627496909 | Lamé 260.0 | thin wall 200.0
b / a 2 | hoop stress at the bore 166.6610191803395 | Lamé 166.66666666666666 | thin wall 100.0
b / a 3 | hoop stress at the bore 124.99506084735357 | Lamé 125.0 | thin wall 50.0
b / a 4 | hoop stress at the bore 113.32847778699471 | Lamé 113.33333333333333 | thin wall 33.333333333333336
b / a 5 | hoop stress at the bore 108.32853722389413 | Lamé 108.33333333333333 | thin wall 25.0

Figure 1: 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:

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))

interference 0.022766078542970976 | contact pressure 10.32905415375535 | hoop stress at both bores 129.48207171314743

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):

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)

shrunk on: contact pressure 10.329562448362418 | formula 10.32905415375535 | tube's bore -37.183777861563506 | formula -37.184594953519266
under pressure: tube's bore 129.47703628309154 | ring's bore 129.4661037940268 | single wall's bore 166.66666666666669

(a) One wall

(b) Tube and shrunk-on ring

Figure 2: 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):

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()

Figure 3: The hoop stress across the wall

See Stress at a shoulder fillet for a stress concentration, and Cantilever deflection for a plane solid against beam theory.