Clamped-clamped beam under a central point load

Statically-indeterminate first-order beam with both ends fully clamped, carrying a single transverse point load \(P\) at midspan. Symmetry plus compatibility at the clamps gives the elementary closed form (Roark & Young 1989, Table 8 case 5; Timoshenko 1955 §40):

\[R_\mathrm{left} = R_\mathrm{right} = \tfrac{P}{2}, \qquad M_\mathrm{end} = \pm \tfrac{P L}{8},\]

with the mid-span deflection (the indeterminate clamping moments stiffen the beam by a factor of 4 vs. the simply-supported case)

\[\delta_\mathrm{mid} = \frac{P\, L^{3}}{192\, E I},\]

and the mid-span bending moment \(M_\mathrm{mid} = +PL/8\) (numerically equal to the end-moment magnitude because the moment diagram is anti-symmetric about midspan).

Implementation

A 20-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam. Both end nodes are clamped (all 6 DOFs); a single \(-P\,\hat y\) force is applied at the midspan node.

The Hermite cubic basis recovers Euler-Bernoulli kinematics exactly under nodal point loads, so every node hits the analytical value to machine precision and the clamping reactions / moments come out cleanly.

References

• Timoshenko, S. P. (1955) Strength of Materials, Part I, 3rd ed., Van Nostrand, §40 (statically indeterminate beams).

• Roark, R. J. and Young, W. C. (1989) Roark’s Formulas for Stress and Strain, 6th ed., McGraw-Hill, Table 8 case 5 (clamped-clamped beam, central point load).

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §2.5.

Vendor cross-references

Source	Reported δ_mid [m]	Problem ID / location
Closed form (Euler-Bernoulli)	5.000 × 10⁻⁵	Timoshenko SoM Part I §5.7
Gere & Goodno (2018) Table 10-1	5.000 × 10⁻⁵	fixed-fixed beam, central load
MAPDL Verification Manual	5.00 × 10⁻⁵	VM-2 Beam stresses and deflections (clamped variant)
Abaqus Verification Manual	5.00 × 10⁻⁵	AVM 1.5.x fixed-fixed beam family

from __future__ import annotations

import numpy as np
import pyvista as pv
from vtkmodules.util.vtkConstants import VTK_LINE

import femorph_solver
from femorph_solver import ELEMENTS

Problem data

E = 2.0e11  # Pa (steel)
NU = 0.30
RHO = 7850.0  # kg/m^3
b = 0.050  # square cross-section side [m]
A = b * b
I = b**4 / 12.0  # noqa: E741
J = 2.0 * I
L = 1.0  # span [m]
P = 5.0e3  # central point load magnitude [N], applied downward (-y)

N_ELEM = 20  # even so a node lands exactly at midspan

# Closed-form quantities -----------------------------------------------

R_published = P / 2.0
M_end_published = P * L / 8.0
delta_mid_published = P * L**3 / (192.0 * E * I)

print(f"P = {P:.1f} N, L = {L:.2f} m, EI = {E * I:.3e} N m^2")
print(f"R_L = R_R = P/2          = {R_published:.4f} N")
print(f"|M_end|  = P L / 8       = {M_end_published:.4f} N m  (clamping moment)")
print(f"δ_mid  = P L^3/(192 EI)  = {delta_mid_published:.4e} m")
print("stiffness ratio C-C / S-S = 4x   (192 / 48 = 4)")

P = 5000.0 N, L = 1.00 m, EI = 1.042e+05 N m^2
R_L = R_R = P/2          = 2500.0000 N
|M_end|  = P L / 8       = 625.0000 N m  (clamping moment)
δ_mid  = P L^3/(192 EI)  = 2.5000e-04 m
stiffness ratio C-C / S-S = 4x   (192 / 48 = 4)

Build a 20-element BEAM2 line

points = np.array(
    [[i * L / N_ELEM, 0.0, 0.0] for i in range(N_ELEM + 1)],
    dtype=np.float64,
)
cells_list: list[int] = []
for i in range(N_ELEM):
    cells_list.extend([2, i, i + 1])
cells = np.asarray(cells_list, dtype=np.int64)
cell_types = np.full(N_ELEM, VTK_LINE, dtype=np.uint8)
grid = pv.UnstructuredGrid(cells, cell_types, points)

m = femorph_solver.Model.from_grid(grid)
m.assign(
    ELEMENTS.BEAM2,
    material={"EX": E, "PRXY": NU, "DENS": RHO},
    real=(A, I, I, J),
)

Boundary conditions: both ends clamped

right = N_ELEM + 1
m.fix(nodes=[1], dof="ALL")
m.fix(nodes=[right], dof="ALL")

Apply the central point load

mid_node = N_ELEM // 2 + 1
m.apply_force(mid_node, fy=-P)

Static solve + reaction extraction

res = m.solve_static()
dof = m.dof_map()


def _react(node: int, dof_idx: int) -> float:
    rows = np.where((dof[:, 0] == node) & (dof[:, 1] == dof_idx))[0]
    return float(res.reaction[rows[0]]) if len(rows) else 0.0


R_left = _react(1, 1)
R_right = _react(right, 1)
M_left = _react(1, 5)  # ROTZ reaction = clamp moment
M_right = _react(right, 5)

print()
print("reactions  (femorph-solver) → (analytical)")
print(f"  R_left_UY    = {R_left:+.4f}  →  {+R_published:+.4f} N")
print(f"  R_right_UY   = {R_right:+.4f}  →  {+R_published:+.4f} N")
print(f"  M_left_RZ    = {M_left:+.4f}  →  {+M_end_published:+.4f} N m")
print(f"  M_right_RZ   = {M_right:+.4f}  →  {-M_end_published:+.4f} N m")

# Vertical equilibrium and exact reaction match
assert np.isclose(R_left + R_right, P, rtol=1e-12), "vertical equilibrium failed"
assert np.isclose(R_left, R_published, rtol=1e-12), "left reaction off"
assert np.isclose(R_right, R_published, rtol=1e-12), "right reaction off"
# Clamp moment-reactions are equal in magnitude, opposite in sign — both
# act to resist the downward load at midspan.  M_left = +PL/8 (CCW about
# +z) lifts the right end against the load; M_right = -PL/8 (CW)
# mirrors it on the other side.
assert np.isclose(M_left, +M_end_published, rtol=1e-12), "left clamp moment off"
assert np.isclose(M_right, -M_end_published, rtol=1e-12), "right clamp moment off"
assert np.isclose(M_left + M_right, 0.0, atol=1e-9), "moment anti-symmetry violated"

reactions  (femorph-solver) → (analytical)
  R_left_UY    = +2500.0000  →  +2500.0000 N
  R_right_UY   = +2500.0000  →  +2500.0000 N
  M_left_RZ    = +625.0000  →  +625.0000 N m
  M_right_RZ   = -625.0000  →  -625.0000 N m

Mid-span deflection

mid_uy = float(res.displacement[np.where((dof[:, 0] == mid_node) & (dof[:, 1] == 1))[0][0]])
err_mid = (mid_uy - (-delta_mid_published)) / (-delta_mid_published)
print()
print(f"δ_mid  computed   = {mid_uy:+.4e} m")
print(f"δ_mid  published  = {-delta_mid_published:+.4e} m")
print(f"relative error    = {err_mid * 100:+.6f} %")
assert abs(err_mid) < 1e-8, "mid-span deflection error too large for Bernoulli BEAM2"

δ_mid  computed   = -2.5000e-04 m
δ_mid  published  = -2.5000e-04 m
relative error    = -0.000000 %

Render the deflected line

pts = np.asarray(grid.points)
disp_y = np.zeros(N_ELEM + 1)
for i in range(N_ELEM + 1):
    rows = np.where((dof[:, 0] == i + 1) & (dof[:, 1] == 1))[0]
    if len(rows):
        disp_y[i] = float(res.displacement[rows[0]])
warped = grid.copy()
warped.points = pts + np.column_stack([np.zeros(N_ELEM + 1), 200.0 * disp_y, np.zeros(N_ELEM + 1)])
warped["uy"] = disp_y

plotter = pv.Plotter(off_screen=True, window_size=(720, 280))
plotter.add_mesh(grid, color="grey", line_width=2, opacity=0.5)
plotter.add_mesh(warped, scalars="uy", line_width=5, cmap="viridis")
plotter.add_points(
    np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]]),
    render_points_as_spheres=True,
    point_size=18,
    color="black",
    label="clamps",
)
plotter.add_points(
    np.array([[0.5 * L, 0.0, 0.0]]),
    render_points_as_spheres=True,
    point_size=14,
    color="#d62728",
    label=f"P = {P:.0f} N",
)
plotter.view_xy()
plotter.camera.zoom(1.05)
plotter.show()