8.14  Beam theory

The assumptions

Four assumptions define the Euler-Bernoulli beam. The deformations are small, so that the linear strain measure 8.7.3 applies. The material is linearly elastic. Poisson effects are neglected, so that the axial stress follows the one-dimensional Hooke’s law \(\sigma = E\varepsilon\) rather than the full 8.8.7. Finally, and this is the assumption that does the work, every cross section remains plane and perpendicular to the neutral axis during deformation.

The last assumption is a statement about geometry, and it is what Figure 8.14.1 illustrates. Longitudinal lines in the beam become curved, while lines that were vertical stay straight and merely rotate. Somewhere between the top and the bottom is a surface that neither stretches nor shortens: the neutral plane, whose intersection with the plane of bending is the neutral axis. For a homogeneous beam the neutral axis passes through the centroid of the cross section, which is why the second moment of area is always taken about the centroid.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 3.6))

L_fig, h = 10.0, 2.0
xs = np.linspace(0, L_fig, 200)

# Undeformed beam with a grid of material lines
ax1.add_patch(plt.Rectangle((0, -h/2), L_fig, h, fill=False, edgecolor='k', lw=1.5))
for yy in np.linspace(-h/2, h/2, 5):
    ax1.plot([0, L_fig], [yy, yy], 'k-', lw=0.4)
for xx in np.linspace(0, L_fig, 11):
    ax1.plot([xx, xx], [-h/2, h/2], 'k-', lw=0.4)
ax1.plot([0, L_fig], [0, 0], 'r--', lw=1.5)
ax1.text(L_fig*0.5, h/2 + 0.4, "Before deformation", ha='center')
ax1.text(L_fig + 0.3, 0, r"neutral axis", color='r', va='center', fontsize=9)

# Deformed beam: neutral axis y(X1), cross sections normal to it
kappa = 0.012
y_of = lambda X: -kappa*X**2
dy_of = lambda X: -2*kappa*X
for xx in np.linspace(0, L_fig, 11):
    th = dy_of(xx)
    nx, ny = -np.sin(th), np.cos(th)     # normal to the deformed axis
    ax2.plot([xx - nx*h/2, xx + nx*h/2],
             [y_of(xx) - ny*h/2, y_of(xx) + ny*h/2], 'k-', lw=0.4)
for s in np.linspace(-h/2, h/2, 5):
    th = dy_of(xs)
    ax2.plot(xs - np.sin(th)*s, y_of(xs) + np.cos(th)*s, 'k-', lw=0.4)
ax2.plot(xs, y_of(xs), 'r--', lw=1.5)
ax2.text(L_fig*0.5, h/2 + 0.4, "After deformation", ha='center')

for a in (ax1, ax2):
    a.set_aspect('equal')
    a.set_xlim(-0.8, L_fig + 2.2)
    a.set_ylim(-h*1.9, h*1.5)
    a.axis('off')
plt.tight_layout()
plt.show()

Kinematics

Let \(X_1\) run along the beam and \(X_2\) across it, with \(X_3\) out of the plane, and let \(y(X_1)\) denote the deflection of the neutral axis. Consider a material point \(P\) that sits at \((X_1^P, X_2^P, X_3^P)\) before deformation. Its cross section translates by \(y\) and rotates by an angle \(\theta\), so after deformation the point sits at

\[ \bm x = \begin{bmatrix} X_1^P - X_2^P \sin\theta\\ y + X_2^P \cos\theta\\ X_3^P \end{bmatrix} \tag{8.14.1}\]

The rotation angle is not independent. Since the cross section stays perpendicular to the neutral axis, it rotates by exactly the slope of that axis, so \(\theta\) is fixed by \(y\) through

\[ \tan\theta = \frac{dy}{dX_1} = y' \]

and for small deformations the small angle approximations \(\sin\theta \approx \theta \approx \tan\theta\) and \(\cos\theta \approx 1\) turn 8.14.1 into

\[ \bm u = \bm x - \bm X = \begin{bmatrix} -X_2^P\,\theta\\ y\\ 0 \end{bmatrix} = \begin{bmatrix} -X_2 \dfrac{dy}{dX_1}\\[8pt] y\\ 0 \end{bmatrix} \tag{8.14.2}\]

which is the complete displacement field of the beam. A three-dimensional vector field has been reduced to one scalar function \(y(X_1)\) of one variable. Everything else in this chapter is bookkeeping applied to 8.14.2.

Strain and stress

The strain follows from the continuum definition 8.7.3 applied to 8.14.2. We let SymPy form the displacement gradient and its symmetric part rather than differentiating nine entries by hand.

X_1, X_2, X_3 = sp.symbols('X_1 X_2 X_3', real=True)
E_mod = sp.symbols('E', positive=True)

y = sp.Function('y')(X_1)

uu = sp.Matrix([-X_2*sp.diff(y, X_1), y, 0])
XX = sp.Matrix([X_1, X_2, X_3])

grad_u = uu.jacobian(XX)
eps = sp.simplify((grad_u + grad_u.T)/2)

ltx(r"\nabla \bm u =", grad_u, r",\qquad \bm\varepsilon = \tfrac12(\nabla\bm u + \nabla\bm u^\mathsf{T}) =", eps)

The displacement gradient has two non-zero entries, but its symmetric part has only one. The off-diagonal terms \(-y'\) and \(+y'\) cancel exactly, which is the algebraic signature of the perpendicularity assumption: the cross section rotates rigidly with the axis, and a rigid rotation produces no strain. What survives is a single axial strain

\[ \boxed{\varepsilon_{11} = -X_2 \frac{d^2 y}{dX_1^2}} \tag{8.14.3}\]

proportional to the distance from the neutral axis and to the curvature of the beam. Fibres above the neutral axis (\(X_2 > 0\)) shorten when the curvature is positive, fibres below it stretch, and the fibre on the axis does neither, which is what earned it the name.

The stress follows from Hooke’s law, and since we have neglected Poisson effects it is the one-dimensional form,

sigma = sp.simplify(E_mod*eps)
ltx(r"\bm\sigma = E\bm\varepsilon =", sigma)

so that the only non-zero stress component is

\[ \boxed{\sigma_{11} = -E X_2 \frac{d^2 y}{dX_1^2}} \tag{8.14.4}\]

The stress varies linearly through the depth of the section, from maximum compression at one surface through zero at the neutral axis to maximum tension at the other. Material near the neutral axis is barely worked at all, which is why I-beams and box sections put their material at the extremes and why a hollow tube is so much more efficient in bending than a solid rod of the same mass.

Internal forces and the second moment of area

The stress distribution 8.14.4 must be statically equivalent to the internal forces on the cut. Its resultant moment about the \(X_3\) axis is obtained by integrating the stress times its lever arm over the cross section,

\[ M = \iint -X_2\,\sigma_{11}\,dX_2\,dX_3 = \iint E X_2^2 \frac{d^2y}{dX_1^2}\,dX_2\,dX_3 = E\frac{d^2y}{dX_1^2} \underbrace{\iint X_2^2\,dX_2\,dX_3}_{I} \]

where the curvature came out of the integral because it does not vary over the cross section. The remaining integral is a purely geometric property,

\[ \boxed{I = \iint X_2^2\, dX_2\, dX_3} \tag{8.14.5}\]

the second moment of area about the neutral axis, with units of length to the fourth power. This gives the moment-curvature relation

\[ \boxed{M = EI \frac{d^2 y}{dX_1^2}} \tag{8.14.6}\]

in which the product \(EI\) appears as a single quantity, the bending stiffness. Material and geometry enter only through their product, so a stiffer material and a deeper section are interchangeable as far as deflection is concerned.

Eliminating the curvature between 8.14.4 and 8.14.6 gives the result that every engineer knows,

\[ \boxed{\sigma_{11} = -\frac{M X_2}{I}} \tag{8.14.7}\]

and evaluating it at the extreme fibre, at distance \(c\) from the neutral axis, gives the maximum bending stress

\[ \sigma_{\max} = \frac{Mc}{I} = \frac{M}{W}, \qquad W := \frac{I}{c} \tag{8.14.8}\]

where \(W\) is the section modulus, the single number that determines how much moment a section can carry at a given stress. Values of \(I\) and \(W\) for standard sections are tabulated in every handbook, and a CAD system reports them for arbitrary sections directly. For a rectangle of width \(b\) and height \(h\) the integral is elementary.

b, h_sec = sp.symbols('b h', positive=True)

I_rect = sp.integrate(sp.integrate(X_2**2, (X_2, -h_sec/2, h_sec/2)), (X_3, 0, b))
W_rect = sp.simplify(I_rect/(h_sec/2))

ltx(r"I_{\text{rect}} &=", I_rect, r"\\ W_{\text{rect}} &=", W_rect, aligned=True)

The height enters \(I\) cubed. Doubling the height of a beam makes it eight times stiffer and four times stronger, while doubling its width does neither more than doubling. This single fact governs almost every decision about how to orient a structural member.

Equilibrium of a slice

8.14.6 relates the moment to the deflection, but nothing yet relates the moment to the applied load. That link comes from equilibrium, applied not to the whole beam but to a slice of length \(dX_1\) carrying a distributed load \(q\) per unit length. The internal forces on the left face of the slice are a normal force \(N\), a shear force \(V\) and a bending moment \(M\), and on the right face each has changed by its own derivative times \(dX_1\).

fig, ax = plt.subplots(figsize=(6.5, 3.4))

wd, ht = 3.0, 2.0
ax.add_patch(plt.Rectangle((0, 0), wd, ht, fill=False, edgecolor='k', lw=1.5))

# distributed load
for xx in np.linspace(0.15, wd-0.15, 9):
    ax.annotate("", xy=(xx, ht), xytext=(xx, ht+0.55),
                arrowprops=dict(arrowstyle='->', lw=1))
ax.plot([0, wd], [ht+0.55, ht+0.55], 'k-', lw=1)
ax.text(wd/2, ht+0.75, r"$q$", ha='center')

# shear
ax.annotate("", xy=(0, ht/2-0.6), xytext=(0, ht/2+0.6), arrowprops=dict(arrowstyle='->', lw=1.6, color='b'))
ax.text(-0.42, ht/2, r"$V$", color='b', va='center')
ax.annotate("", xy=(wd, ht/2+0.6), xytext=(wd, ht/2-0.6), arrowprops=dict(arrowstyle='->', lw=1.6, color='b'))
ax.text(wd+0.12, ht/2, r"$V + \frac{\partial V}{\partial X_1}dX_1$", color='b', va='center', fontsize=9)

# moment
ax.text(0.28, ht/2+0.02, r"$M$", color='g', va='center')
ax.text(wd-1.05, ht/2+0.02, r"$M + \frac{\partial M}{\partial X_1}dX_1$", color='g', va='center', fontsize=9)

ax.plot([0, wd], [-0.35, -0.35], 'k-', lw=0.8)
ax.plot([0, 0], [-0.25, -0.45], 'k-', lw=0.8)
ax.plot([wd, wd], [-0.25, -0.45], 'k-', lw=0.8)
ax.text(wd/2, -0.75, r"$dX_1$", ha='center')

ax.set_xlim(-1.0, wd+1.6)
ax.set_ylim(-1.1, ht+1.2)
ax.set_aspect('equal')
ax.axis('off')
plt.tight_layout()
plt.show()

Vertical force equilibrium of the slice gives the change in shear force as the applied load, and moment equilibrium about the right face gives the change in bending moment as the shear force,

\[ \frac{dV}{dX_1} = q, \qquad \frac{dM}{dX_1} = V \tag{8.14.9}\]

These two relations are worth committing to memory in their own right, because they say that the shear diagram is the integral of the load and the moment diagram is the integral of the shear. A moment diagram has a maximum wherever the shear crosses zero, and the shear jumps by exactly the size of any point load.

The elastic line equation

Now assemble the pieces. Substituting the moment-curvature relation 8.14.6 into 8.14.9 twice gives

\[ \boxed{EI \frac{d^4 y}{dX_1^4} = q(X_1)} \tag{8.14.10}\]

the elastic line equation, a fourth-order ordinary differential equation for the deflection. The three-dimensional continuum problem has been reduced to this. Everything that made the general problem hard, the coupling between displacement components and the two extra spatial dimensions, was eliminated by the kinematic assumption.

Because 8.14.10 is obtained by differentiating 8.14.6, the intermediate derivatives are all quantities we have names for. Integrating downward from the load gives a ladder,

\[ \begin{array}{llll} EI\,\dfrac{d^4y}{dX_1^4} = q(X_1) & \quad\text{load} \\[8pt] EI\,\dfrac{d^3y}{dX_1^3} = V(X_1) & \quad\text{shear force} \\[8pt] EI\,\dfrac{d^2y}{dX_1^2} = M(X_1) & \quad\text{bending moment} \\[8pt] \dfrac{dy}{dX_1} = \theta(X_1) & \quad\text{slope} \\[8pt] y(X_1) & \quad\text{deflection} \end{array} \tag{8.14.11}\]

so a single solution of 8.14.10 produces the shear diagram, the moment diagram, the rotation and the deflection at once, together with the bending stress through 8.14.7. There is no need to draw shear and moment diagrams by hand as a separate exercise: they are derivatives of the answer.

Boundary conditions

A fourth-order equation needs four conditions, two at each end. At each end there are four quantities available, and exactly two of them may be prescribed: the deflection, the slope, the moment and the shear force. Which two depends on the support.

\[ \begin{array}{ll} y(0) = 0 & \qquad y(L) = 0\\[4pt] \dfrac{dy}{dX_1}(0) = 0 & \qquad \dfrac{dy}{dX_1}(L) = 0\\[8pt] \dfrac{d^2y}{dX_1^2}(0) = \dfrac{M}{EI} & \qquad \dfrac{d^2y}{dX_1^2}(L) = \dfrac{M}{EI}\\[8pt] \dfrac{d^3y}{dX_1^3}(0) = \dfrac{V}{EI} & \qquad \dfrac{d^3y}{dX_1^3}(L) = \dfrac{V}{EI} \end{array} \tag{8.14.12}\]

The pairing follows a simple rule: a support that prevents a motion prescribes the corresponding kinematic quantity, and one that permits it prescribes the corresponding force quantity instead. A clamped end prevents both translation and rotation, so \(y = 0\) and \(y' = 0\). A pinned or roller support prevents translation but permits rotation, so \(y = 0\) and \(M = 0\). A free end prevents nothing, so \(M\) and \(V\) are prescribed, taking the values of whatever external moment and force are applied there, and zero if none are.

fig, axes = plt.subplots(1, 3, figsize=(11, 2.8))

def beam_line(a):
    a.plot([0, 4], [0, 0], 'k-', lw=3)
    a.set_xlim(-1.2, 4.6); a.set_ylim(-1.5, 1.3)
    a.set_aspect('equal'); a.axis('off')

# clamped
beam_line(axes[0])
axes[0].add_patch(plt.Rectangle((-0.35, -0.9), 0.35, 1.8, color='k'))
for yy in np.linspace(-0.85, 0.85, 7):
    axes[0].plot([-0.7, -0.35], [yy-0.2, yy], 'k-', lw=0.8)
axes[0].set_title(r"Clamped" "\n" r"$y=0,\;y'=0$", fontsize=10)

# pinned
beam_line(axes[1])
axes[1].plot([0, -0.45, 0.45, 0], [0, -0.75, -0.75, 0], 'k-', lw=1.2)
for xx in np.linspace(-0.5, 0.5, 6):
    axes[1].plot([xx, xx-0.2], [-0.75, -0.98], 'k-', lw=0.8)
axes[1].plot([-0.6, 0.6], [-0.75, -0.75], 'k-', lw=1.2)
axes[1].set_title(r"Pinned" "\n" r"$y=0,\;M=0$", fontsize=10)

# free
beam_line(axes[2])
axes[2].annotate("", xy=(4, -0.9), xytext=(4, 0.2), arrowprops=dict(arrowstyle='->', lw=1.5, color='b'))
axes[2].text(4.1, -0.5, r"$V$", color='b')
axes[2].set_title(r"Free" "\n" r"$M$ and $V$ prescribed", fontsize=10)

plt.tight_layout()
plt.show()

Example 1: A cantilever under its own distributed load

A cantilever of length \(L\) is clamped at \(X_1 = 0\) and free at \(X_1 = L\), carrying a uniformly distributed downward load of magnitude \(q_0\). Determine the deflection, and the tip deflection in particular.

The clamped end contributes \(y(0) = 0\) and \(y'(0) = 0\). The free end carries no applied moment and no applied force, so \(M(L) = 0\) and \(V(L) = 0\). With \(q = -q_0\) for a downward load, the problem is fully specified.

x, L, EI = sp.symbols('x L EI', positive=True)
q_0, E_s, I_s = sp.symbols('q_0 E I', positive=True)
w = sp.Function('w')

general = sp.dsolve(sp.Eq(E_s*I_s*sp.diff(w(x), x, 4), -q_0), w(x)).rhs
constants = sorted(general.free_symbols - {x, L, E_s, I_s, q_0}, key=lambda s: s.name)

M_expr = E_s*I_s*sp.diff(general, x, 2)
V_expr = E_s*I_s*sp.diff(general, x, 3)

bcs = [general.subs(x, 0),               # y(0) = 0
       sp.diff(general, x).subs(x, 0),   # y'(0) = 0
       M_expr.subs(x, L),                # M(L) = 0
       V_expr.subs(x, L)]                # V(L) = 0

sol = sp.solve(bcs, constants, dict=True)[0]
y_sol = sp.simplify(general.subs(sol))

ltx(r"y(x) &=", y_sol,
    r"\\ y(L) &=", sp.simplify(y_sol.subs(x, L)), aligned=True)

The tip deflection is \(-q_0 L^4/(8EI)\), negative because the beam sags in the direction the load pushes it, and matching the tabulated elementary case. The fourth power of the length is the reason cantilevers are so sensitive to their span: doubling the reach of an overhang increases its droop sixteenfold.

The whole ladder 8.14.11 follows by differentiation, so the shear and moment diagrams cost nothing extra.

theta_sol = sp.diff(y_sol, x)
M_sol = E_s*I_s*sp.diff(y_sol, x, 2)
V_sol = E_s*I_s*sp.diff(y_sol, x, 3)

ltx(r"\theta(x) &=", sp.simplify(theta_sol),
    r"\\ M(x) &=", sp.simplify(M_sol),
    r"\\ V(x) &=", sp.simplify(V_sol), aligned=True)

Both diagrams are exact and both make sense at a glance. The moment is largest in magnitude at the clamp, where it equals \(-q_0L^2/2\), the moment of the whole load about the support. The shear at the clamp is \(-q_0 L\), the entire applied load, and both fall to zero at the free tip as the boundary conditions demanded.

Piecewise beams and matching conditions

8.14.10 assumes that \(q(X_1)\) is a smooth function of position, and real beams rarely oblige. A point load, an intermediate support, a step in cross section or an applied couple all break the beam into segments, within each of which the load is smooth but across which the solution is not.

The remedy is to write one differential equation per segment and to join the solutions with matching conditions, which the Swedish literature calls skarvvillkor. Each additional segment brings four more constants and therefore needs four more conditions, and the interface supplies exactly four. Two are kinematic: the beam does not break apart and does not develop a kink, so the deflection and the slope are continuous. Two are static: the internal moment and shear are continuous except for any concentrated moment or force applied at the joint, which enters as a jump.

For a joint at \(X_1 = a\) carrying a downward point load \(P\) and no applied couple, the four conditions are

\[ \begin{array}{ll} y_1(a) = y_2(a) & \text{no separation}\\[4pt] y_1'(a) = y_2'(a) & \text{no kink}\\[4pt] M_1(a) = M_2(a) & \text{no applied couple}\\[4pt] V_1(a) - V_2(a) = P & \text{shear jumps by the point load} \end{array} \tag{8.14.13}\]

The last line is the one to check. Applying 8.14.13 to a simply supported beam of span \(L\) with a central point load should reproduce the textbook deflection \(PL^3/(48EI)\).

P = sp.symbols('P', positive=True)
w1f, w2f = sp.Function('w_1'), sp.Function('w_2')

gen1 = sp.dsolve(sp.Eq(E_s*I_s*sp.diff(w1f(x), x, 4), 0), w1f(x)).rhs
gen2 = sp.dsolve(sp.Eq(E_s*I_s*sp.diff(w2f(x), x, 4), 0), w2f(x)).rhs

C_1 = sorted(gen1.free_symbols - {x, E_s, I_s}, key=lambda s: s.name)
gen2 = gen2.subs({c: sp.Symbol('D' + c.name[1:]) for c in C_1})
C_2 = sorted(gen2.free_symbols - {x, E_s, I_s}, key=lambda s: s.name)

M_1, M_2 = E_s*I_s*sp.diff(gen1, x, 2), E_s*I_s*sp.diff(gen2, x, 2)
V_1, V_2 = E_s*I_s*sp.diff(gen1, x, 3), E_s*I_s*sp.diff(gen2, x, 3)

a = L/2
equations = [gen1.subs(x, 0), M_1.subs(x, 0),                        # pinned at 0
             gen2.subs(x, L), M_2.subs(x, L),                        # pinned at L
             (gen1 - gen2).subs(x, a),                               # continuity of y
             (sp.diff(gen1, x) - sp.diff(gen2, x)).subs(x, a),       # continuity of slope
             (M_1 - M_2).subs(x, a),                                 # continuity of moment
             (V_1 - V_2).subs(x, a) - P]                             # shear jump

sol = sp.solve(equations, C_1 + C_2, dict=True)[0]
y_mid = sp.simplify(gen1.subs(sol).subs(x, a))

ltx(r"y\!\left(\tfrac{L}{2}\right) =", y_mid)

The result is \(-PL^3/(48EI)\), downward and of exactly the tabulated magnitude, so the sign convention in 8.14.13 is the right one.

Example 2: A propped cantilever with an interior load

We now solve a case that no elementary table covers. A beam is pinned at \(X_1 = 0\), clamped at the far end, and carries a downward point load \(P\) where its two segments meet, at \(X_1 = L_1\). The second segment has length \(L_2\). The structure is statically indeterminate, so equilibrium alone cannot find the reactions and the deflection must be solved for first.

L_1, L_2 = sp.symbols('L_1 L_2', positive=True)

gen1 = sp.dsolve(sp.Eq(E_s*I_s*sp.diff(w1f(x), x, 4), 0), w1f(x)).rhs
gen2 = sp.dsolve(sp.Eq(E_s*I_s*sp.diff(w2f(x), x, 4), 0), w2f(x)).rhs
C_1 = sorted(gen1.free_symbols - {x, E_s, I_s}, key=lambda s: s.name)
gen2 = gen2.subs({c: sp.Symbol('D' + c.name[1:]) for c in C_1})
C_2 = sorted(gen2.free_symbols - {x, E_s, I_s}, key=lambda s: s.name)

M_1, M_2 = E_s*I_s*sp.diff(gen1, x, 2), E_s*I_s*sp.diff(gen2, x, 2)
V_1, V_2 = E_s*I_s*sp.diff(gen1, x, 3), E_s*I_s*sp.diff(gen2, x, 3)
L_tot = L_1 + L_2

equations = [gen1.subs(x, 0), M_1.subs(x, 0),                        # pinned at x = 0
             gen2.subs(x, L_tot), sp.diff(gen2, x).subs(x, L_tot),   # clamped at the far end
             (gen1 - gen2).subs(x, L_1),
             (sp.diff(gen1, x) - sp.diff(gen2, x)).subs(x, L_1),
             (M_1 - M_2).subs(x, L_1),
             (V_1 - V_2).subs(x, L_1) - P]

sol = sp.solve(equations, C_1 + C_2, dict=True)[0]
y_1 = sp.simplify(gen1.subs(sol))
y_2 = sp.simplify(gen2.subs(sol))

ltx(r"y_1(x) =", y_1)

With the symbolic solution in hand we insert a specific design. Take a square steel section \(100 \times 100~\text{mm}\), so that \(I = bh^3/12\) and the extreme fibre sits at \(c = 50~\text{mm}\), with \(L_1 = 1000~\text{mm}\), \(L_2 = 2000~\text{mm}\) and \(P = 1000~\text{N}\).

values = {E_s: 210e3, I_s: 100**4/12, P: 1000.0, L_1: 1000.0, L_2: 2000.0}
c_fibre = 50.0

pieces = []
for y_seg, (x_a, x_b) in [(y_1, (sp.Integer(0), L_1)), (y_2, (L_1, L_tot))]:
    pieces.append(dict(
        y=sp.lambdify(x, y_seg.subs(values), 'numpy'),
        th=sp.lambdify(x, sp.diff(y_seg, x).subs(values), 'numpy'),
        M=sp.lambdify(x, (E_s*I_s*sp.diff(y_seg, x, 2)).subs(values), 'numpy'),
        V=sp.lambdify(x, (E_s*I_s*sp.diff(y_seg, x, 3)).subs(values), 'numpy'),
        span=(float(x_a.subs(values)), float(x_b.subs(values)))))

def sample(key):
    xs, vs = [], []
    for p in pieces:
        xi = np.linspace(*p['span'], 200)
        vi = np.asarray(p[key](xi), dtype=float) * np.ones_like(xi)
        xs.append(xi); vs.append(vi)
    return np.concatenate(xs), np.concatenate(vs)

x_v, y_v = sample('y')
_,   M_v = sample('M')
_,   V_v = sample('V')
sigma_v = M_v * c_fibre / values[I_s]

ltx(r"y_{\max} &=", float(np.min(y_v)), r"~\text{mm}",
    r"\\ M_{\max} &=", float(np.max(np.abs(M_v))), r"~\text{Nmm}",
    r"\\ \sigma_{\max} &=", float(np.max(np.abs(sigma_v))), r"~\text{MPa}", aligned=True)

The largest deflection is \(0.14~\text{mm}\) downward and the largest bending stress is \(3.1~\text{MPa}\), so this section is enormously oversized for the load, which is what we would expect from a \(100~\text{mm}\) square carrying a kilonewton. Plotting the four diagrams together shows the structure of the solution.

fig, axes = plt.subplots(4, 1, figsize=(8, 9), sharex=True)

for a, v, lbl, col in zip(axes, (y_v, V_v, M_v, sigma_v),
                          (r"$y$ [mm]", r"$V$ [N]", r"$M$ [Nmm]", r"$\sigma$ [MPa]"),
                          ('C0', 'C4', 'C2', 'C3')):
    a.plot(x_v, v, color=col, lw=1.8)
    a.fill_between(x_v, 0, v, color=col, alpha=0.25)
    a.axhline(0, color='k', lw=0.8)
    a.axvline(float(values[L_1]), color='k', ls=':', lw=1)
    a.set_ylabel(lbl)
    a.grid(True, alpha=0.4)

axes[-1].set_xlabel(r"$X_1$ [mm]")
axes[0].set_title("Deflection, shear, moment and bending stress")
plt.tight_layout()
plt.show()

The shear is piecewise constant, as it must be when no distributed load acts, and it jumps by exactly \(P\) at the joint. The moment is piecewise linear, peaks at the joint where the shear changes sign, and is negative at the clamp because the beam curves the other way there. The bending stress is the moment scaled by \(c/I\), so it is the moment diagram in different units, and it passes through zero at the inflection point where the curvature reverses. That point is where the top surface of the beam stops being in compression and starts being in tension, which matters when the two surfaces are made of different materials or when only one of them is welded.

Transverse shear stress

The Euler-Bernoulli kinematics produced a strain tensor with a single non-zero entry, so it predicts \(\varepsilon_{12} = 0\) and hence, through \(\tau = G\gamma\), a transverse shear stress of exactly zero. Yet 8.14.9 plainly contains a shear force \(V\), and a force cannot be transmitted by a stress that is zero. The two statements contradict each other.

The contradiction is a consequence of the perpendicularity assumption, which forbids the cross section from shearing while the equilibrium of the beam requires it to. The standard resolution is to accept the kinematics for computing the bending stress and to obtain the shear stress from equilibrium instead. Take a slice of length \(\Delta X_1\) and cut it horizontally at some height, keeping the portion above the cut. The axial stress 8.14.7 acts on both ends of that portion, and because the moment changes across the slice the two are not equal. The imbalance must be carried by a shear stress \(\sigma_{21}\) on the horizontal cut.

Writing the two axial resultants and the shear force on the horizontal face,

\[ F_1 = \iint -\frac{M + \Delta M}{I} X_2\, dX_2\, dX_3, \qquad F_2 = \iint \frac{M}{I} X_2\, dX_2\, dX_3, \qquad F_3 = \sigma_{21}\,\Delta X_1\, b \]

where \(b\) is the width of the section at the height of the cut, equilibrium in the axial direction gives \(F_1 + F_2 + F_3 = 0\), which reduces to

\[ \iint -\frac{\Delta M}{I} X_2\, dX_2\, dX_3 - \sigma_{21}\Delta X_1 b = 0 \]

Using \(\Delta M \approx (dM/dX_1)\Delta X_1 = V \Delta X_1\) from 8.14.9 and cancelling \(\Delta X_1\) leaves

\[ \boxed{\tau = \sigma_{12} = \sigma_{21} = -\frac{V Q}{I b}}, \qquad Q := \iint X_2\, dX_2\, dX_3 \tag{8.14.14}\]

where \(Q\) is the first moment of the area lying beyond the cut, taken about the neutral axis. The symmetry \(\sigma_{12} = \sigma_{21}\) that we proved in Chapter 8.7 is what allows the horizontal shear computed here to be read as the vertical shear on the cross section.

The distribution is instructive. At the top and bottom surfaces the area beyond the cut vanishes, so \(Q = 0\) and the shear stress is zero, as it must be on a free surface. At the neutral axis \(Q\) is largest, so the shear stress peaks there, exactly where the bending stress vanishes. For a rectangular section the profile is parabolic with a maximum of \(\tfrac{3}{2}V/A\), fifty percent above the average.

xi, V_sym = sp.symbols('xi V', real=True)

# first moment of the area above the cut at height xi, about the neutral axis
Q_rect = b * sp.integrate(X_2, (X_2, xi, h_sec/2))
tau_rect = sp.simplify(V_sym*Q_rect/(I_rect*b))

ltx(r"Q(\xi) &=", sp.simplify(Q_rect),
    r"\\ \tau(\xi) &=", tau_rect,
    r"\\ \tau(0) &=", sp.simplify(tau_rect.subs(xi, 0)), aligned=True)

With \(A = bh\) the peak value is \(\tfrac{3}{2}V/A\) as claimed. In a slender beam the bending stress is much larger than this, which is why shear rarely governs the design of a long beam and often governs a short deep one.

The Timoshenko-Ehrenfest beam

The contradiction in the previous section is resolved properly by weakening the kinematic assumption. Keep the requirement that cross sections remain plane, but drop the requirement that they remain perpendicular to the neutral axis. The rotation of the cross section then becomes an independent function \(\psi\), related to the slope of the axis by

\[ \psi = y' + \varphi \]

where \(\varphi\) is the additional shear rotation. The displacement field 8.14.2 becomes

\[ \bm u = \begin{bmatrix} -X_2\psi \\ y \\ 0 \end{bmatrix} \tag{8.14.15}\]

and the same continuum machinery applied to it now yields a strain tensor with a shear entry that no longer cancels.

psi = sp.Function('psi')(X_1)
y_t = sp.Function('y')(X_1)

uu_t = sp.Matrix([-X_2*psi, y_t, 0])
grad_t = uu_t.jacobian(XX)
eps_t = sp.simplify((grad_t + grad_t.T)/2)

ltx(r"\nabla\bm u =", grad_t, r",\qquad \bm\varepsilon =", eps_t)

The axial strain is now \(\varepsilon_{11} = -X_2\,d\psi/dX_1\), driven by the gradient of the section rotation rather than by the curvature of the axis, and the shear strain is

\[ 2\varepsilon_{12} = \frac{dy}{dX_1} - \psi = -\varphi =: -\gamma \]

Integrating the stresses over the cross section as before gives the two internal force relations

\[ M = EI\frac{d\psi}{dX_1}, \qquad V = kAG\gamma \tag{8.14.16}\]

in which \(k\) is a shear correction factor. It is needed because the assumption of a plane cross section forces the shear stress to be uniform through the depth, whereas 8.14.14 says it is parabolic. The factor \(k\) smears the real distribution into an equivalent uniform one and takes values around \(5/6\) for a rectangle. It is not an exact quantity, and different authors derive different values from different equivalence criteria.

Substituting 8.14.16 into the slice equilibrium 8.14.9 gives the governing pair

\[ EI\frac{d^3\psi}{dX_1^3} = q, \qquad \frac{dy}{dX_1} = \psi - \frac{EI}{kAG}\frac{d^2\psi}{dX_1^2} \tag{8.14.17}\]

The first equation is solved for \(\psi\) and the second then integrated for \(y\), which introduces a fourth constant, so the Timoshenko beam still needs four boundary conditions. The deflection comes out as the Euler-Bernoulli answer plus a shear correction proportional to \(EI/(kAG)\).

That group is the whole story. Write it as

\[ \frac{EI}{kAG L^2} \sim \frac{E}{G}\cdot\frac{I}{A L^2} \sim \frac{E}{G}\left(\frac{r}{L}\right)^2 \]

where \(r = \sqrt{I/A}\) is the radius of gyration of the section, a length of the order of its depth. Since \(E/G = 2(1+\nu) \approx 2.6\) by 8.8.5, the correction is a few times the square of the depth-to-span ratio. A beam ten times longer than it is deep gets a correction of a few percent and the Euler-Bernoulli model is entirely adequate. A beam three times longer than it is deep gets a correction approaching thirty percent, and a beam as deep as it is long is not a beam at all and should be modelled as a two-dimensional continuum.

When the beam model applies

Every result in this chapter descends from 8.14.2, so every result inherits its limitations. The model assumes that the cross section moves rigidly, which fails near a point where load is introduced, near a support, near a hole and near an abrupt change of section, because in those regions the cross section warps. Saint-Venant’s principle says that the disturbance decays within a distance of the order of the section depth, so the beam solution is accurate everywhere except in short zones around the discontinuities. Those zones are exactly where a finite element model earns its cost, and they are where the stress concentrations of a real design live.

The model also assumes bending in one plane about a principal axis of the section. Loading a section obliquely, or loading an unsymmetric section at all, couples bending about the two axes and, for open thin-walled sections, couples bending to torsion as well.

What we have gained is worth the restrictions. A partial differential equation for a vector field over a volume became an ordinary differential equation for a scalar on a line, solvable in closed form, differentiable into every quantity a designer needs, and cheap enough to put inside an optimisation loop. That trade, a kinematic assumption in exchange for a dimensional reduction, is the pattern behind every structural theory, and recognising it as a special case of the continuum model of Chapter 8.8 rather than as a separate subject is what makes the pattern visible.