---
jupyter:
  description: 'The catenary arch with planex2d and NGSolve: a catenary and a semicircular arch of the same span and
    rise under their own weight; the catenary stays in compression, the semicircle opens at its haunches and its crown,
    and the catenary''s thrust tends to w a as the arch gets thinner.'
---

# 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.

![The largest principal stress in the two arches under their own weight, tension in red](../_static/gallery/catenary-arch.png){.result}

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

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

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

```python
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()})
```

```python nbsphinx="hide-input"
Row(*[(face, f"The {name}", {"dashed": [centerlines[name]]}) for name, face in arches.items()],
    caption="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.

```python
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:

```python
for name, (mesh, s1, forces) in results.items():
    print(name, "| supports", forces["vertical"], "N, weight", rho * g * arches[name].mass, "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:

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

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:

```python
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")
```

```python nbsphinx="hidden"
def sampled(mesh, cf):
    """sigma 1 in kPa for Row, at points inside the arch"""
    def field(face, plot=True, **_):
        (u0, v0), (u1, v1) = face.bounding_box
        gu, gv = np.meshgrid(np.linspace(u0, u1, 301), np.linspace(v0, v1, 151))
        uv = np.c_[gu.ravel(), gv.ravel()]
        uv = uv[face.Distance(uv, signed=True) < -0.005]
        values = np.asarray(cf(mesh(uv[:, 0], uv[:, 1]))).ravel() / 1e3
        return {"points": [tuple(p) for p in uv], "values": values, "label": "σ1, kPa"}
    return field
```

```python nbsphinx="hide-input"
Row(*[(arches[name], f"The {name}", {"field": sampled(results[name][0], results[name][1]), "cmap": MAP.reversed(),
                                     "clim": (-60, 60)}) for name in arches],
    caption="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.

```python nbsphinx="hidden"
import sys
import matplotlib.pyplot as plt
sys.path.insert(0, "..")
from gallery import save_figure
fig, axes = plt.subplots(1, 2, figsize=(9, 3.6))
for ax, name in zip(axes, ("catenary", "semicircle")):
    mesh, s1, forces = results[name]
    px.PlotField(s1 / 1e3, mesh, arches[name], ax=ax, colormap=MAP.reversed(), min=-60, max=60, samples=600)
save_figure(fig, "catenary-arch", bare=True, panels=[0, 1], bars={3: "σ1, kPa"})
plt.close(fig)
```