The catenary arch

A chain hanging between two points takes the shape of a catenary, \(y = a \cosh(x/a)\), and is in pure tension. Turned upside down, the same shape carries its own weight in pure compression, which is all that masonry can carry. Robert Hooke’s rule says it in one line, “as hangs the flexible line, so but inverted will stand the rigid arch”. Other shapes need their thickness to keep the line of thrust inside them, and where it leaves, they open.

The example compares a catenary and a semicircular arch of the same span and rise, both 0.5 m thick, under their own weight in 2D elasticity. The supports must carry the weight exactly. The catenary’s horizontal thrust must tend to \(w \cdot a\), its weight per length times its parameter, as the arch gets thinner. The catenary must stay in compression where the semicircle opens.

import math
import numpy as np
from scipy.optimize import brentq
import ngsolve as ngs
import planex as px

Two arches

Both arches span 10 m and rise 5 m. The semicircle has the radius 5, and the catenary \(y = 5 - a\,(\cosh(x/a) - 1)\) reaches the ground at \(\pm 5\) for the parameter \(a\) that solves \(5 = a\,(\cosh(5/a) - 1)\). The catenary’s centerline is drawn with px.Curve from this formula, with the curve’s parameter \(t\) running along \(x\) and the numbers written into the formula by an f-string. The semicircle’s centerline is an arc. Each centerline, offset by half the thickness to both sides, is a strip, cut off at the ground by a large rectangle whose bottom line is named support (Figure 1):

span, rise, t = 10.0, 5.0, 0.5
a = brentq(lambda a: a * (math.cosh(span / 2 / a) - 1) - rise, 1, 20)
centerlines = {
    "catenary": px.Curve(("t", f"{rise} - {a} * (cosh(t / {a}) - 1)"), (-span / 2, span / 2)),
    "semicircle": px.QuickSketch().Move((span / 2, 0), dir=(0, 1)).Arc(span / 2, 180).Wires()[0],
}
ground = px.QuickSketch((-20, 0)).Line(40, name="support").Line(0, 20).Line(-40, 0).Close().Face()


def arch(centerline, thickness):
    return (centerline.Offset(thickness / 2) * ground).Face()


arches = {name: arch(line, t) for name, line in centerlines.items()}
print("a", a, "| areas", {name: face.mass for name, face in arches.items()})

a 3.0937961388712787 | areas {'catenary': 7.477808753322162, 'semicircle': 7.853981633974485}

(a) The catenary

(b) The semicircle

Figure 1: A catenary and a semicircular arch of 10 m span and 5 m rise, 0.5 m thick

Their own weight

The arches are concrete, 2400 kg/m³ with \(E = 20\,\text{GPa}\) and \(\nu = 0.2\), in plane strain, as a long vault is, so the stresses are those of a slice 1 m deep. The supports hold the arches’ feet. solve returns the displacement, the stress and the largest principal stress \(\sigma_1\), which is positive where the material is in tension. The supports’ reactions come from the solved system itself. Its residual at the supports’ degrees of freedom, paired with a unit displacement along \(x\) at the left support, is the force the left support takes along \(x\), and the same along \(y\) at both supports is the vertical force they take.

rho, g, E, nu = 2400.0, 9.81, 20e9, 0.2
mu, lam = E / 2 / (1 + nu), E * nu / (1 + nu) / (1 - 2 * nu)


def solve(face, maxh):
    mesh = ngs.Mesh(px.OCCGeometry(face).GenerateMesh(maxh=maxh))
    mesh.Curve(3)
    V = ngs.VectorH1(mesh, order=3, dirichlet="support")
    u, v = V.TnT()

    def stress(q):
        eps = ngs.Sym(ngs.Grad(q))
        return 2 * mu * eps + lam * ngs.Trace(eps) * ngs.Id(2)

    stiffness = ngs.BilinearForm(ngs.InnerProduct(stress(u), ngs.Sym(ngs.Grad(v))) * ngs.dx).Assemble()
    weight = ngs.LinearForm(-rho * g * v[1] * ngs.dx).Assemble()
    displacement = ngs.GridFunction(V)
    displacement.vec.data = stiffness.mat.Inverse(V.FreeDofs(), inverse="umfpack") * weight.vec
    s = stress(displacement)
    s1 = (s[0, 0] + s[1, 1]) / 2 + ngs.sqrt(((s[0, 0] - s[1, 1]) / 2) ** 2 + s[0, 1] ** 2)
    residual = displacement.vec.CreateVector()
    residual.data = stiffness.mat * displacement.vec - weight.vec

    def reaction(direction):
        unit = ngs.GridFunction(V)
        unit.Set(direction, definedon=mesh.Boundaries("support"))
        return ngs.InnerProduct(residual, unit.vec)

    forces = {"thrust": reaction(ngs.CF((ngs.IfPos(-ngs.x, 1, 0), 0))), "vertical": reaction(ngs.CF((0, 1)))}
    return mesh, s1, forces


results = {name: solve(face, t / 5) for name, face in arches.items()}

The supports must carry the arches’ weight, \(\rho g\) times their area:

for name, (mesh, s1, forces) in results.items():
    print(name, "| supports", forces["vertical"], "N, weight", rho * g * arches[name].mass, "N")

catenary | supports 176057.56885455034 N, weight 176057.529288217 N
semicircle | supports 184914.14353504227 N, weight 184914.14359029525 N

The thrust

An arch pushes its supports apart, and the horizontal force at a support is its thrust. In a thin catenary arch, it is the same everywhere along the arch, \(H = w \cdot a\), with \(w = \rho g t\) its weight per length. A real arch has a thickness, and the thinner it is, the closer it comes to the thin-arch value. The catenary’s thrust over \(w \cdot a\) is printed for thicknesses of 0.5, 0.25 and 0.125:

for thickness in (0.5, 0.25, 0.125):
    mesh, s1, forces = solve(arch(centerlines["catenary"], thickness), thickness / 5)
    w = rho * g * thickness
    print("thickness", thickness, "| thrust", forces["thrust"], "N, w a", w * a, "N, ratio", forces["thrust"] / (w * a))

thickness 0.5 | thrust 35967.0851653051 N, w a 36420.16814679269 N, ratio 0.9875595582189125
thickness 0.25 | thrust 18153.13716039782 N, w a 18210.084073396345 N, ratio 0.9968727814342317
thickness 0.125 | thrust 9097.92255948727 N, w a 9105.042036698173 N, ratio 0.9992180731091402

Each halving of the thickness divides the difference by four, so the thin-arch value is the limit, approached with the square of the thickness.

Where the arches open

The largest principal stress \(\sigma_1\) shows where an arch is in tension (Figure 2). Sampled just inside the outer and the inner face, by the angle from the right support, the semicircle is in tension on the outside at its haunches and on the inside at its crown. For the catenary, the cell prints its largest \(\sigma_1\) above the supports’ first meter, against its compression:

mesh, s1, forces = results["semicircle"]
for face, radius in (("outer", span / 2 + t / 2 - 0.01), ("inner", span / 2 - t / 2 + 0.01)):
    values = {angle: s1(mesh(radius * math.cos(math.radians(angle)), radius * math.sin(math.radians(angle)))) / 1e3
              for angle in range(10, 91, 5)}
    tension = [angle for angle, value in values.items() if value > 1]
    peak = max(values, key=values.get)
    print("semicircle,", face, "face: in tension from", min(tension), "to", max(tension), "degrees, at most",
          round(values[peak]), "kPa at", peak, "degrees")
mesh, s1, forces = results["catenary"]
above = [s1(mesh(*v.point)) / 1e3 for v in mesh.vertices if v.point[1] > 1]
print("catenary above 1 m: sigma 1 from", round(min(above), 1), "to", round(max(above), 2), "kPa")

semicircle, outer face: in tension from 20 to 50 degrees, at most 181 kPa at 35 degrees
semicircle, inner face: in tension from 70 to 90 degrees, at most 224 kPa at 90 degrees
catenary above 1 m: sigma 1 from -0.2 to 0.24 kPa

(a) The catenary

(b) The semicircle

Figure 2: The largest principal stress under the arches' own weight, tension in red (red from 60 kPa on)

The semicircle opens where its line of thrust leaves it, outside at the haunches and inside at the crown, where a masonry arch of this shape cracks, and inside at its feet, where the supports hold it. The catenary has no such places. Its largest principal stress stays at about zero across its thickness, the smaller one carries the compression along the arch, and the whole arch is in compression.