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.

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

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:

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:

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

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:

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

fixed and free | critical load 345.41 | Euler 345.44
pinned and pinned | critical load 1379.93 | Euler 1381.74
fixed and pinned | critical load 2819.2 | Euler 2827.96
fixed and fixed | critical load 5500.7 | Euler 5526.98

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:

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

fixed and free | with shear 345.35
pinned and pinned | with shear 1380.33
fixed and pinned | with shear 2822.03
fixed and fixed | with shear 5504.38

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

(a) Fixed and free (rotated 90 degrees)

(b) Pinned and pinned (rotated 90 degrees)

(c) Fixed and pinned (rotated 90 degrees)

(d) Fixed and fixed (rotated 90 degrees)

Figure 1: The buckled shapes, the displacement drawn at a tenth of the length