Note
Go to the end to download the full example code.
Cantilever with off-tip point load — load-position shape function#
A clamped-free prismatic beam of length \(L\) with a single transverse point load \(P\) applied at an interior position \(x = a \le L\). The Euler-Bernoulli closed form (Timoshenko 1955 §32; Roark Table 8 case 1) splits at the load:
so the deflection grows as \(x^{2}\,(3 a - x)\) while the load is “above” \(x\) and as \(a^{2}\,(3 x - a)\) past the load (a straight ramp of slope \(P\,a^{2} / (2\, E\, I)\) away from the load). Three closed-form check points fall out:
Deflection at the load
\[v(a) = \frac{P\, a^{3}}{3\, E\, I}.\]Tip deflection \(v(L)\)
\[v(L) = \frac{P\, a^{2}\,(3 L - a)}{6\, E\, I}.\]Tip slope \(\theta(L)\) (= the slope past the load)
\[\theta(L) = \frac{P\, a^{2}}{2\, E\, I}.\]
For the tip-load limit \(a = L\) the formula collapses to \(v(L) = P L^{3} / (3 E I)\), recovering the Cantilever tip deflection — Euler-Bernoulli closed form result. For \(a = L/2\) (midspan loading) the tip deflection is \(5 P L^{3} / (48 E I) = 5/16\) of the tip-loaded value — load position matters far more than load magnitude on a slender cantilever.
Implementation#
A 40-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam. The load is applied at the closest mesh node to \(x = a\); the example reports the FE deflection / slope at three check points along with the closed-form expectation for each.
Hermite-cubic shape functions interpolate Bernoulli kinematics exactly on the polynomial moment field, so the FE response under this point load matches the closed form to machine precision at the loaded node and the tip — well within a 0.1 % tolerance.
References#
Timoshenko, S. P. (1955) Strength of Materials, Part I: Elementary Theory and Problems, 3rd ed., Van Nostrand, §32 (cantilever with intermediate load).
Roark, R. J. and Young, W. C. (1989) Roark’s Formulas for Stress and Strain, 6th ed., McGraw-Hill, Table 8 case 1 (cantilever, concentrated load, intermediate position).
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 (Hermite-cubic shape functions).
from __future__ import annotations
import numpy as np
import pyvista as pv
import femorph_solver
from femorph_solver import ELEMENTS
Problem data — slender beam, midspan load#
E = 2.0e11
NU = 0.30
RHO = 7850.0
WIDTH = 0.05
HEIGHT = 0.05
A_section = WIDTH * HEIGHT
I_z = WIDTH * HEIGHT**3 / 12.0
I_y = HEIGHT * WIDTH**3 / 12.0
J = (1.0 / 3.0) * min(WIDTH, HEIGHT) ** 3 * max(WIDTH, HEIGHT)
L = 1.0 # span [m]
a = L / 2.0 # load position from the clamp [m]
P = 1.0e3 # transverse load [N, downward = -y]
EI = E * I_z
def v_closed_form(x: float) -> float:
"""Closed-form deflection v(x) for the off-tip cantilever load."""
if x <= a:
return P * x**2 * (3.0 * a - x) / (6.0 * EI)
return P * a**2 * (3.0 * x - a) / (6.0 * EI)
def theta_closed_form(x: float) -> float:
"""Closed-form slope dv/dx(x) for the off-tip cantilever load."""
if x <= a:
return P * x * (2.0 * a - x) / (2.0 * EI)
return P * a**2 / (2.0 * EI)
v_at_load_pub = v_closed_form(a)
v_tip_pub = v_closed_form(L)
theta_tip_pub = theta_closed_form(L)
print("Cantilever with off-tip point load")
print(f" L = {L} m, load at x = a = {a} m (a / L = {a / L:.3f})")
print(f" E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2, P = {P} N")
print()
print("Closed-form references (Roark Table 8 case 1):")
print(f" v(a) = P a^3 / (3 E I) = {v_at_load_pub * 1e3:+.4e} mm")
print(f" v(L) = P a^2 (3 L - a) / (6 E I) = {v_tip_pub * 1e3:+.4e} mm")
print(f" theta(L) = P a^2 / (2 E I) = {theta_tip_pub:+.4e} rad")
Cantilever with off-tip point load
L = 1.0 m, load at x = a = 0.5 m (a / L = 0.500)
E = 2.00e+11 Pa, I = 5.208e-07 m^4, EI = 1.042e+05 N m^2, P = 1000.0 N
Closed-form references (Roark Table 8 case 1):
v(a) = P a^3 / (3 E I) = +4.0000e-01 mm
v(L) = P a^2 (3 L - a) / (6 E I) = +1.0000e+00 mm
theta(L) = P a^2 / (2 E I) = +1.2000e-03 rad
Build a 40-element BEAM2 line#
N_ELEM = 40
xs = np.linspace(0.0, L, N_ELEM + 1)
pts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))
cells = np.empty((N_ELEM, 3), dtype=np.int64)
cells[:, 0] = 2
cells[:, 1] = np.arange(N_ELEM)
cells[:, 2] = np.arange(1, N_ELEM + 1)
grid = pv.UnstructuredGrid(
cells.ravel(),
np.full(N_ELEM, 3, dtype=np.int64),
pts,
)
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.BEAM2,
material={"EX": E, "PRXY": NU, "DENS": RHO},
real=(A_section, I_z, I_y, J),
)
Boundary conditions: full clamp at x = 0#
m.fix(nodes=1, dof="ALL")
Apply the load at the closest node to x = a#
Suppress out-of-plane motion at every node so the strictly-2D response matches the closed-form Bernoulli kinematics.
Solve + extract#
res = m.solve_static()
dof_map = m.dof_map()
disp = np.asarray(res.displacement, dtype=np.float64)
def _node_dof(node_id: int, dof_idx: int) -> float:
"""0=UX, 1=UY, 2=UZ, 3=ROTX, 4=ROTY, 5=ROTZ."""
rows = np.where(dof_map[:, 0] == node_id)[0]
for r in rows:
if int(dof_map[r, 1]) == dof_idx:
return float(disp[r])
return 0.0
v_at_load_fe = _node_dof(i_load + 1, 1)
v_tip_fe = _node_dof(N_ELEM + 1, 1)
theta_tip_fe = _node_dof(N_ELEM + 1, 5) # ROTZ — in-plane rotation
# The applied load points in -y (gravity-like), so FE deflections
# come out negative. The closed-form formulas above return the
# magnitude. Compare on absolute value across the board so the
# sign convention drops out.
err_load = (abs(v_at_load_fe) - v_at_load_pub) / v_at_load_pub * 100.0
err_tip = (abs(v_tip_fe) - v_tip_pub) / v_tip_pub * 100.0
err_theta = (abs(theta_tip_fe) - theta_tip_pub) / theta_tip_pub * 100.0
print()
print(f"{'quantity':<22} {'FE':>14} {'published':>14} {'err':>9}")
print(f"{'-' * 22:<22} {'-' * 14:>14} {'-' * 14:>14} {'-' * 9:>9}")
print(
f"{'v(a)':<22} {v_at_load_fe * 1e3:>10.4f} mm "
f"{v_at_load_pub * 1e3:>10.4f} mm {err_load:>+8.4f}%"
)
print(f"{'v(L)':<22} {v_tip_fe * 1e3:>10.4f} mm {v_tip_pub * 1e3:>10.4f} mm {err_tip:>+8.4f}%")
print(
f"{'|theta(L)|':<22} {abs(theta_tip_fe):>11.6e} {theta_tip_pub:>11.6e} {err_theta:>+8.4f}%"
)
assert abs(err_load) < 0.1, f"v(a) deviation {err_load:.4f} % exceeds 0.1 %"
assert abs(err_tip) < 0.1, f"v(L) deviation {err_tip:.4f} % exceeds 0.1 %"
assert abs(err_theta) < 0.1, f"theta(L) deviation {err_theta:.4f} % exceeds 0.1 %"
quantity FE published err
---------------------- -------------- -------------- ---------
v(a) -0.4000 mm 0.4000 mm -0.0000%
v(L) -1.0000 mm 1.0000 mm -0.0000%
|theta(L)| 1.200000e-03 1.200000e-03 -0.0000%
Plot the deflection profile against the closed form#
fe_uy = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)])
analytical_uy = np.array([v_closed_form(x) for x in xs])
import matplotlib.pyplot as plt # noqa: E402
fig, ax = plt.subplots(1, 1, figsize=(7.0, 4.0))
xs_dense = np.linspace(0.0, L, 401)
ax.plot(
xs_dense,
[v_closed_form(x) * 1e3 for x in xs_dense],
"-",
color="#1f77b4",
lw=2.0,
label="closed form (Roark Table 8 case 1)",
)
ax.plot(xs, fe_uy * 1e3, "o", color="#d62728", markersize=5, label="BEAM2 FE")
ax.axvline(a, color="grey", ls="--", lw=1.0, label=f"load at x = a = {a:.3f} m")
ax.set_xlabel("x [m]")
ax.set_ylabel("v(x) [mm] (downward deflection negative)")
ax.set_title("Cantilever with off-tip point load — FE vs closed form", fontsize=11)
ax.legend(loc="lower left", fontsize=9)
ax.grid(True, ls=":", alpha=0.5)
fig.tight_layout()
fig.show()
Render the deformed shape#
g = m.grid.copy()
disp_xyz = np.zeros((g.n_points, 3))
for ni in range(g.n_points):
disp_xyz[ni, 1] = _node_dof(int(ni + 1), 1)
g.point_data["displacement"] = disp_xyz
g.point_data["UY"] = disp_xyz[:, 1]
scale = 1.0 / max(abs(v_tip_fe), 1e-12) * 0.05
warped = g.copy()
warped.points = np.asarray(g.points) + scale * disp_xyz
plotter = pv.Plotter(off_screen=True, window_size=(900, 320))
plotter.add_mesh(g, color="grey", opacity=0.4, line_width=2, label="undeformed")
plotter.add_mesh(warped, scalars="UY", cmap="coolwarm", line_width=4)
plotter.add_points(
np.array([[0.0, 0.0, 0.0]]),
render_points_as_spheres=True,
point_size=18,
color="black",
label="clamp",
)
plotter.add_points(
np.array([[a, scale * v_at_load_fe, 0.0]]),
render_points_as_spheres=True,
point_size=14,
color="#d62728",
label=f"load (P = {P:.0f} N at x = {a:.2f} m)",
)
plotter.view_xy()
plotter.add_legend()
plotter.show()
Take-aways#
print()
print("Take-aways:")
print(f" • v(a) within {abs(err_load):.4f} % of P a³ / (3 E I).")
print(f" • v(L) within {abs(err_tip):.4f} % of P a² (3 L − a) / (6 E I).")
print(f" • |theta(L)| within {abs(err_theta):.4f} % of P a² / (2 E I).")
ratio = v_tip_pub / (P * L**3 / (3.0 * EI))
print(
f" • Compared to a tip-loaded cantilever of the same beam, the midspan-loaded "
f"tip deflects only {ratio:.4f} of the tip-loaded value (= 5/16 for a = L/2)."
)
Take-aways:
• v(a) within 0.0000 % of P a³ / (3 E I).
• v(L) within 0.0000 % of P a² (3 L − a) / (6 E I).
• |theta(L)| within 0.0000 % of P a² / (2 E I).
• Compared to a tip-loaded cantilever of the same beam, the midspan-loaded tip deflects only 0.3125 of the tip-loaded value (= 5/16 for a = L/2).
Total running time of the script: (0 minutes 0.283 seconds)