"""Simply-supported beam under central point load.
Reference geometry
------------------
Prismatic slender beam of length :math:`L`, rectangular cross-
section :math:`w \\times h`, isotropic elastic material. Simply
supported at both ends — pinned at :math:`x = 0` (all three
translations fixed), roller at :math:`x = L` (transverse
translations fixed, axial free) — and loaded by a central
transverse point load :math:`P` at :math:`x = L/2`.
Closed-form quantities
----------------------
For a slender Euler-Bernoulli beam, integrating
:math:`EI\\,w'''' = 0` piecewise over the two spans
:math:`[0, L/2]` and :math:`[L/2, L]` with the central point-load
jump in the shear ``V`` gives
.. math::
\\delta_{\\text{mid}} = \\frac{P L^{3}}{48 E I}, \\qquad
I = \\frac{w h^{3}}{12}.
The support reactions are each :math:`R = P/2` by symmetry.
References
----------
* Timoshenko, S. P. *Strength of Materials, Part I*, 3rd ed.,
Van Nostrand, 1955, §5.6 — simply-supported beam under a
concentrated mid-span load.
* Gere & Goodno (2018), *Mechanics of Materials* 9th ed., §9.3
Table 9-2, case 5.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.5.x** simply-supported-beam family.
* **NAFEMS** *Background to Benchmarks* (§3.2, simply-supported
beam with central load).
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any
import numpy as np
import pyvista as pv
import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.validation._benchmark import (
BenchmarkProblem,
PublishedValue,
)
[docs]
@dataclass
class SimplySupportedBeamCentralLoad(BenchmarkProblem):
"""Simply-supported beam mid-span deflection under central point load."""
name: str = "ss_beam_central_load"
description: str = (
"Prismatic simply-supported beam with central transverse "
"point load — Euler-Bernoulli closed-form "
"(Timoshenko 1955 §5.6)."
)
analysis: str = "static"
#: Beam length [m].
L: float = 1.0
#: Cross-section width [m].
width: float = 0.05
#: Cross-section height [m].
height: float = 0.05
#: Young's modulus [Pa].
E: float = 2.0e11
#: Poisson's ratio.
nu: float = 0.3
#: Density [kg/m^3] — irrelevant for the static solution.
rho: float = 7850.0
#: Total transverse mid-span load [N] (acts in ``-z``).
P: float = 1.0e3
@property
def published_values(self) -> tuple[PublishedValue, ...]:
I = self.width * self.height**3 / 12.0 # noqa: E741
delta = self.P * self.L**3 / (48.0 * self.E * I)
return (
PublishedValue(
name="mid_span_deflection",
value=delta,
unit="m",
source="Timoshenko 1955 §5.6",
formula="delta = P L^3 / (48 E I)",
tolerance=0.05,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
# Requires an even ``nx`` so a single node sits exactly at
# the mid-span — this is where the concentrated load lands.
nx = int(mesh_params.get("nx", 40))
if nx % 2 != 0:
nx += 1
ny = int(mesh_params.get("ny", 3))
nz = int(mesh_params.get("nz", 3))
xs = np.linspace(0.0, self.L, nx + 1)
ys = np.linspace(0.0, self.width, ny + 1)
zs = np.linspace(0.0, self.height, nz + 1)
grid = pv.StructuredGrid(
*np.meshgrid(xs, ys, zs, indexing="ij")
).cast_to_unstructured_grid()
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.HEX8(integration="enhanced_strain"),
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
)
pts = np.asarray(m.grid.points)
# Simply-supported BCs on a 3D hex beam require knife-edge
# support lines (not point supports) so the problem doesn't
# degenerate into a locally-loaded 3D stress concentration
# at a single node. Both support lines sit at the bottom
# (z = 0) face, running the full width in y.
#
# Left end (x = 0, z = 0, all y): pin UZ along the line
# + pin UX and UY at one
# single node to suppress
# the six rigid-body DOFs.
# Right end (x = L, z = 0, all y): pin UZ along the line
# (roller — UX free).
left_line = np.where((pts[:, 0] < 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in left_line:
m.fix(nodes=int(n + 1), dof="UZ", value=0.0)
# Single-node axial + transverse pin at the (0, 0, 0) corner
# to remove rigid-body X and Y modes without over-constraining.
pin_node = np.where((pts[:, 0] < 1e-9) & (pts[:, 1] < 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in pin_node:
m.fix(nodes=int(n + 1), dof="UX", value=0.0)
m.fix(nodes=int(n + 1), dof="UY", value=0.0)
# An extra UY pin at the (L, 0, 0) corner kills the last
# rigid-body Y-translation mode without blocking the axial
# (UX) strain that the roller support must allow.
right_pin_node = np.where(
(pts[:, 0] > self.L - 1e-9) & (pts[:, 1] < 1e-9) & (pts[:, 2] < 1e-9)
)[0]
for n in right_pin_node:
m.fix(nodes=int(n + 1), dof="UY", value=0.0)
# Right-end roller — UZ fixed along the full-width line.
right_line = np.where((pts[:, 0] > self.L - 1e-9) & (pts[:, 2] < 1e-9))[0]
for n in right_line:
m.fix(nodes=int(n + 1), dof="UZ", value=0.0)
# Central point load: lump onto the mid-span bottom-line
# nodes so the support line and load line share the same
# cross-section level, giving the cleanest discrete mirror
# of the Euler-Bernoulli problem.
x_mid = 0.5 * self.L
load_nodes = np.where((np.abs(pts[:, 0] - x_mid) < 1e-9) & (pts[:, 2] < 1e-9))[0]
fz_per_node = -self.P / len(load_nodes)
for n in load_nodes:
m.apply_force(int(n + 1), fz=fz_per_node)
# Nodes at the mid-span top face (x = L/2, z = h). The
# top face sits opposite the concentrated-load line on
# the bottom; its deflection is the clean Euler-Bernoulli
# beam-theory deflection without the local 3D stress-
# concentration indentation that contaminates the loaded
# face on a solid mesh.
m._bench_midspan_top_nodes = np.where(
(np.abs(pts[:, 0] - x_mid) < 1e-9) & (pts[:, 2] > self.height - 1e-9)
)[0] # type: ignore[attr-defined]
return m