Simply-supported beam under a central point load#

Classical statically-determinate first-order beam: simple supports (pin + roller) at \(x = 0\) and \(x = L\), with a single transverse point load \(P\) at midspan. Symmetry plus equilibrium gives the elementary closed form

\[R_\mathrm{left} = R_\mathrm{right} = \tfrac{P}{2}, \qquad M_\mathrm{max} = \tfrac{P L}{4}\;\text{at}\;x = L/2,\]

with mid-span and quarter-span deflections (Roark & Young 1989, Table 8 case 1; Timoshenko 1955 §32)

\[\delta(L/2) = \frac{P L^{3}}{48\, E I}, \qquad \delta(L/4) = \frac{11\, P L^{3}}{768\, E I}.\]

Implementation#

A 20-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam. Boundary conditions:

node 1 (pin): UY = UZ = UX = 0; ROTX, ROTY suppressed; ROTZ free (slope rotates).
node \(N_e + 1\) (roller): UY = UZ = 0; ROTX, ROTY suppressed; UX and ROTZ free.
tip load \(-P\,\hat y\) at the central node \(N_e/2 + 1\).

The Hermite-cubic basis recovers Euler-Bernoulli kinematics exactly under nodal point loads, so the FE solution at every node equals the analytical value to machine precision.

References#

Timoshenko, S. P. (1955) Strength of Materials, Part I, 3rd ed., Van Nostrand, §32 (concentrated transverse loads).
Roark, R. J. and Young, W. C. (1989) Roark’s Formulas for Stress and Strain, 6th ed., McGraw-Hill, Table 8 case 1 (simply-supported 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.
Gere, J. M. and Goodno, B. J. (2018) Mechanics of Materials, 9th ed., Cengage, §9.3 Table 9-2 case 5.

Vendor cross-references#

Source	Reported δ_mid [m]	Problem ID / location
Closed form (Euler-Bernoulli)	2.000 × 10⁻⁴	Timoshenko SoM Part I §5.6
Gere & Goodno (2018) Table 9-2	2.000 × 10⁻⁴	SS beam with concentrated mid-load
MAPDL Verification Manual	2.00 × 10⁻⁴	VM-2 Beam stresses and deflections (SS)
Abaqus Verification Manual	2.00 × 10⁻⁴	AVM 1.5.x SS beam family
NAFEMS Background to Benchmarks	2.00 × 10⁻⁴	§3.2 SS beam with central load

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  # tip load magnitude [N], applied downward (-y)

N_ELEM = 20  # must be even so a node lands at midspan

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

R_left_published = P / 2.0
R_right_published = P / 2.0
delta_mid_published = P * L**3 / (48.0 * E * I)
delta_quarter_published = 11.0 * P * L**3 / (768.0 * E * I)
M_max_published = P * L / 4.0

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_left_published:.4f} N")
print(f"M_max  = P L / 4         = {M_max_published:.4f} N m at midspan")
print(f"δ(L/2) = P L^3/(48 EI)   = {delta_mid_published:.4e} m")
print(f"δ(L/4) = 11 P L^3/(768EI) = {delta_quarter_published:.4e} m")

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_max  = P L / 4         = 1250.0000 N m at midspan
δ(L/2) = P L^3/(48 EI)   = 1.0000e-03 m
δ(L/4) = 11 P L^3/(768EI) = 6.8750e-04 m

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#

Standard simply-supported beam:

left pin — UY, UZ, UX = 0; ROTX, ROTY suppressed; ROTZ free.
right roller — UY, UZ = 0; ROTX, ROTY suppressed; UX & ROTZ free.

# Left pin — pin axial too so the system is statically determinate.
m.fix(nodes=[1], dof="UX")
m.fix(nodes=[1], dof="UY")
m.fix(nodes=[1], dof="UZ")
m.fix(nodes=[1], dof="ROTX")
m.fix(nodes=[1], dof="ROTY")

# Right roller — only UY pinned (vertical reaction); UX free.
right = N_ELEM + 1
m.fix(nodes=[right], dof="UY")
m.fix(nodes=[right], dof="UZ")
m.fix(nodes=[right], dof="ROTX")
m.fix(nodes=[right], dof="ROTY")

Apply the central point load#

A single \(-P\hat y\) force at the midspan node. No distributed-load consistent vector needed — the Hermite cubic basis captures the response of a point load on a Bernoulli beam exactly at every node.

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)

print()
print("reactions  (femorph-solver) → (analytical)")
print(f"  R_left_UY   = {R_left:+.4f}  →  {+R_left_published:+.4f} N")
print(f"  R_right_UY  = {R_right:+.4f}  →  {+R_right_published:+.4f} N")

# Vertical equilibrium and exact agreement.
assert np.isclose(R_left + R_right, P, rtol=1e-12), "vertical equilibrium failed"
assert np.isclose(R_left, R_left_published, rtol=1e-12), "left reaction off"
assert np.isclose(R_right, R_right_published, rtol=1e-12), "right reaction off"

reactions  (femorph-solver) → (analytical)
  R_left_UY   = +2500.0000  →  +2500.0000 N
  R_right_UY  = +2500.0000  →  +2500.0000 N

Mid-span and quarter-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"δ(L/2) computed   = {mid_uy:+.4e} m")
print(f"δ(L/2) 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"

quarter_node = N_ELEM // 4 + 1
quarter_uy = float(res.displacement[np.where((dof[:, 0] == quarter_node) & (dof[:, 1] == 1))[0][0]])
err_quarter = (quarter_uy - (-delta_quarter_published)) / (-delta_quarter_published)
print()
print(f"δ(L/4) computed   = {quarter_uy:+.4e} m")
print(f"δ(L/4) published  = {-delta_quarter_published:+.4e} m")
print(f"relative error    = {err_quarter * 100:+.6f} %")
assert abs(err_quarter) < 1e-8, "quarter-span deflection error too large for Bernoulli BEAM2"

δ(L/2) computed   = -1.0000e-03 m
δ(L/2) published  = -1.0000e-03 m
relative error    = -0.000000 %

δ(L/4) computed   = -6.8750e-04 m
δ(L/4) published  = -6.8750e-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",
)
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()

Total running time of the script: (0 minutes 0.171 seconds)

Gallery generated by Sphinx-Gallery