---
jupyter:
  description: 'Euler buckling with planex2d and NGSolve: a slender steel column as a linear buckling eigenvalue problem
    for the four classic end conditions, against Euler''s critical load and Engesser''s shear correction.'
---

# Euler buckling

At what load does a slender column buckle? Euler's answer is the critical load

$$P_{cr} = \frac{\pi^2 E I}{(K L)^2},$$

with $E$ the modulus, $I$ the second moment of area of the section and $K L$ the effective length: $K = 2$ for a
column fixed at its foot and free at its top, 1 for one pinned at both ends, 0.699 for one fixed at the foot and
pinned at the top, and 0.5 for one fixed at both ends. Below this load the straight column is stable, and at it a bent
shape carries the load as well.

Here, NGSolve finds the critical load of a steel strip in plane stress as an eigenvalue problem, for all four end
conditions, and checks it against Euler's formula.

![The buckled shapes of a column for the four classic end conditions](../_static/gallery/euler-buckling.png){.result}

```python
import math

import ngsolve as ngs
import numpy as np
import planex as px
import scipy.sparse as sp
from scipy.sparse.linalg import eigsh
```

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

## The column

The column is a strip of steel 1000 long and 20 wide, with $E = 210000$ in N/mm² and Poisson's ratio 0.3, so it is 50
times as long as it is wide. Its section, per millimeter of depth, has $I = 20^3/12$. The ends are named for the
end conditions:

```python
L, h, E, nu = 1000.0, 20.0, 210e3, 0.3
strip = px.QuickSketch((0, L / 2)).RectangleC(h, L).Face()
strip.edges(near=(0, 0)).name = "bottom"
strip.edges(near=(0, L)).name = "top"
mesh = ngs.Mesh(px.OCCGeometry(strip).GenerateMesh(maxh=h / 4))
```

## Buckling as an eigenvalue problem

Under the axial force $P$, the column carries the stress $\sigma_{yy} = -P/A$. A bent shape $u$ adds the elastic
stiffness $a(u, v)$ and the stiffness of the stress, $g(u, v) = \int \nabla u\, \sigma_0 : \nabla v$, which is negative
in compression. The column buckles at the smallest load factor $\lambda$ for which $(a + \lambda g)(u, v) = 0$ has a
solution, an eigenvalue problem. `g` is assembled for a unit force:

```python
mu, lam = E / (2 * (1 + nu)), E * nu / (1 - nu**2)
V = ngs.VectorH1(mesh, order=3)
u, v = V.TnT()


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


unit = ngs.CF((0, 0, 0, -1 / h), dims=(2, 2))
a = ngs.BilinearForm(ngs.InnerProduct(stress(u), ngs.Sym(ngs.Grad(v))) * ngs.dx).Assemble()
g = ngs.BilinearForm(ngs.InnerProduct(ngs.Grad(u) * unit, ngs.Grad(v)) * ngs.dx).Assemble()
A = sp.csr_matrix(a.mat.CSR(), shape=(V.ndof, V.ndof))
G = sp.csr_matrix(g.mat.CSR(), shape=(V.ndof, V.ndof))
```

## The end conditions

An end condition holds some of the unknowns at zero. A fixed end holds all of an end's edge. A pinned end holds the
middle of the end only, both ways at the foot and sideways at the top, so the end can turn. The unknowns of $x$ come
first in `VectorH1`, then those of $y$:

```python
def edge(name):
    return sorted({d for el in V.Elements(ngs.BND) if el.mat == name for d in el.dofs})


def middle(y):
    k = min(range(mesh.nv), key=lambda i: math.dist(mesh.vertices[i].point, (0, y)))
    return list(V.GetDofNrs(ngs.NodeId(ngs.VERTEX, k)))


sideways = V.ndof // 2
ends = {"fixed and free": (2, edge("bottom")),
        "pinned and pinned": (1, middle(0) + [d for d in middle(L) if d < sideways]),
        "fixed and pinned": (0.699, edge("bottom") + [d for d in middle(L) if d < sideways]),
        "fixed and fixed": (0.5, edge("bottom") + edge("top"))}
```

## The critical loads

For each end condition, SciPy's `eigsh` solves the problem on the unknowns left free. Inverted as
$-g\,x = \mu\,a\,x$, the critical load is $1/\mu$ for the largest $\mu$. The loads are per millimeter of depth, in N:

```python
I = h**3 / 12
modes = {}
for name, (K, held) in ends.items():
    free = np.ones(V.ndof, dtype=bool)
    free[held] = False
    mus, vectors = eigsh(-G[free][:, free], k=1, M=A[free][:, free], which="LA")
    modes[name] = np.zeros(V.ndof)
    modes[name][free] = vectors[:, 0]
    print(name, "| critical load", round(1 / mus[0], 2), "| Euler", round(math.pi**2 * E * I / (K * L) ** 2, 2))
```

The loads lie within half a percent of Euler's, and the shorter the effective length, the lower they lie. The rest is
the shear deformation that Euler's beam leaves out. Engesser's correction,
$P = P_E / (1 + P_E / (\kappa G A))$ with $\kappa = 5/6$ for a rectangle, brings all four within a tenth of a
percent:

```python
for name, (K, held) in ends.items():
    euler = math.pi**2 * E * I / (K * L) ** 2
    print(name, "| with shear", round(euler / (1 + euler / (5 / 6 * mu * h)), 2))
```

Figure 1 shows the four buckled shapes, the sideways displacement drawn at a tenth of the length, turned to lie flat.

```python nbsphinx="hide-input"
mode = ngs.GridFunction(V)
ys = np.linspace(0, L, 201)


def buckled(name):
    mode.vec.FV().NumPy()[:] = modes[name]
    shift = np.array([mode(mesh(0.0, y))[0] for y in ys])
    shift *= 0.1 * L / shift[np.argmax(np.abs(shift))]
    outline = [(-h / 2 + s, y) for s, y in zip(shift, ys)] + [(h / 2 + s, y) for s, y in zip(shift[::-1], ys[::-1])]
    return px.Face([px.FromPolyline(outline, closed=True)])


Row(*[(buckled(name), name.capitalize(), {"rotate": 90}) for name in ends],
    caption="The buckled shapes, the displacement drawn at a tenth of the length")
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
import matplotlib.pyplot as plt
from docs_figures import FACE, INK
from gallery import save_figure
fig, axes = plt.subplots(1, 4, figsize=(6, 4.2))
for ax, name in zip(axes, ends):
    buckled(name).Plot(ax=ax, face_col=FACE, line_col=INK, legend=False)
    ax.set_aspect("equal", adjustable="box")
    ax.set_xlim(-0.15 * L, 0.15 * L)
    ax.axis("off")
save_figure(fig, "euler-buckling", bare=True)
```