Cantilever deflection

How short can a beam be before beam theory fails? Euler–Bernoulli beam theory gives the tip deflection of a cantilever under an end load \(F\) as

\[\delta = \frac{F L^3}{3 E I},\]

assuming that every section stays plane and square to the bent axis, so that the beam does not shear. Timoshenko’s beam theory lets it shear and adds \(F L / (\kappa G A)\), with \(\kappa = 5/6\) for a rectangular section. For a long beam the shear term is small, but for a short one it is not.

The example models the beam as what it is, a plane solid, for lengths from 1 to 20 times its height. NGSolve solves plane stress on a planex2d face, and the tip deflection is compared with both theories and the bending stress away from the ends with \(M c / I\).

This page is part of a collection of physics examples using planex2d.

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

The beam

The beam is a steel strip 20 high, with \(E\) = 210000 N/mm² and Poisson’s ratio 0.3, per millimeter of width. It is held along its whole left end. The load \(F\) = 1000 N acts on its right end, spread over the height as the parabolic shear that beam theory assumes, \(\tau = \frac{3F}{2h}\left(1 - \left(\frac{2y - h}{h}\right)^2\right)\). The function returns the downward deflection of the middle of the end, the bending stress at the top halfway along, and the displacement:

E, nu, h, F = 210e3, 0.3, 20, 1000
mu, lam = E / (2 * (1 + nu)), E * nu / (1 - nu**2)


def cantilever(L):
    """the tip deflection, the bending stress at the top at L/2, and the displacement of a cantilever L long"""
    strip = px.QuickSketch().Rectangle(L, h).Face()
    strip.edges.Min(px.x).name = "clamp"
    strip.edges.Max(px.x).name = "end"
    mesh = ngs.Mesh(px.OCCGeometry(strip).GenerateMesh(maxh=h / 4))
    V = ngs.VectorH1(mesh, order=4, dirichlet="clamp")
    u, v = V.TnT()

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

    a = ngs.BilinearForm(ngs.InnerProduct(stress(u), ngs.Sym(ngs.Grad(v))) * ngs.dx).Assemble()
    shear = -1.5 * F / h * (1 - ((2 * ngs.y - h) / h) ** 2)
    f = ngs.LinearForm(shear * v[1] * ngs.ds("end")).Assemble()
    w = ngs.GridFunction(V)
    w.vec.data = a.mat.Inverse(V.FreeDofs(), inverse="sparsecholesky") * f.vec
    return -w(mesh(L, h / 2))[1], stress(w)[0, 0](mesh(L / 2, h)), (strip, mesh, w)

The deflection

The section’s second moment of area is \(I = h^3/12\), its area \(A = h\) and the shear modulus \(G = E / (2(1 + \nu))\). The plane solid’s deflection is compared with the two theories for lengths from 1 to 20 times the height (Figure 1):

I, A, G = h**3 / 12, h, mu
lengths = [1, 1.5, 2, 3, 4, 6, 8, 12, 20]
tips, tops, beams = {}, {}, {}
for n in lengths:
    L = n * h
    tips[n], tops[n], beams[n] = cantilever(L)
    euler = F * L**3 / (3 * E * I)
    timoshenko = euler + F * L / (5 / 6 * G * A)
    print("L / h", n, "| to Euler-Bernoulli", round(tips[n] / euler, 5), "| to Timoshenko", round(tips[n] / timoshenko, 5))
# plt.plot(lengths, [tips[n] / (F * (n * h)**3 / (3 * E * I)) for n in lengths]); plt.show()

L / h 1 | to Euler-Bernoulli 1.73828 | to Timoshenko 0.97656
L / h 1.5 | to Euler-Bernoulli 1.32592 | to Timoshenko 0.9846
L / h 2 | to Euler-Bernoulli 1.18161 | to Timoshenko 0.9888
L / h 3 | to Euler-Bernoulli 1.07899 | to Timoshenko 0.99293
L / h 4 | to Euler-Bernoulli 1.0434 | to Timoshenko 0.9949
L / h 6 | to Euler-Bernoulli 1.01835 | to Timoshenko 0.99675
L / h 8 | to Euler-Bernoulli 1.00978 | to Timoshenko 0.99762
L / h 12 | to Euler-Bernoulli 1.00387 | to Timoshenko 0.99846
L / h 20 | to Euler-Bernoulli 1.00104 | to Timoshenko 0.9991

Figure 1: The tip deflection against beam theory

At 20 times its height, the beam bends as Euler–Bernoulli says, to a tenth of a percent. At four times, shear adds 4 %, and a beam as long as it is high bends 74 % more. Timoshenko’s shear term covers almost all of it, within 2.3 % even for the shortest beam. The rest comes from the clamp, which holds the whole end face flat where the theory lets the section warp, so the short beams come out a little stiffer.

The stress away from the ends

Away from the clamp and the load, the bending stress follows the beam formula \(\sigma = M c / I\), with the moment \(M = F L / 2\) halfway along and \(c = h / 2\), as Saint-Venant’s principle says. Only in the shortest beams, where the middle is close to both ends, does it differ:

for n in lengths:
    print("L / h", n, "| stress at the top at L/2", round(tops[n], 4), "| M c / I", F * n * h / 2 * (h / 2) / I)
# px.Plot(beams[8][0])

L / h 1 | stress at the top at L/2 149.2181 | M c / I 150.0
L / h 1.5 | stress at the top at L/2 224.6842 | M c / I 225.0
L / h 2 | stress at the top at L/2 299.9604 | M c / I 300.0
L / h 3 | stress at the top at L/2 450.0019 | M c / I 450.0
L / h 4 | stress at the top at L/2 600.0001 | M c / I 600.0
L / h 6 | stress at the top at L/2 900.0 | M c / I 900.0
L / h 8 | stress at the top at L/2 1200.0 | M c / I 1200.0
L / h 12 | stress at the top at L/2 1800.0 | M c / I 1800.0
L / h 20 | stress at the top at L/2 3000.0 | M c / I 3000.0

(a) Twice as long as high

(b) Eight times as long as high

Figure 2: Two cantilevers under their end load, the motion magnified, the unloaded beam dashed

See netgen and NGSolve for meshing a planex2d face, and Euler buckling for the same strip under an axial load.