r""" Propped cantilever with central point load ========================================== A clamped-free cantilever of length :math:`L` propped on a simple support at its tip and loaded by a single concentrated transverse force :math:`P` at midspan :math:`x = L/2`. Statically indeterminate to degree one — the redundant reaction is the upward prop force at the tip. Solving by the force method (or reading Roark Table 8 case 13a) gives the four checkpoint quantities: * **Prop reaction** (upward at :math:`x = L`): .. math:: R_{\mathrm{prop}} = \tfrac{5}{16}\, P. * **Root reaction** (upward at :math:`x = 0`): .. math:: R_{\mathrm{root}} = \tfrac{11}{16}\, P. * **Root bending moment** (negative, hogging): .. math:: M_{\mathrm{root}} = -\tfrac{3}{16}\, P\, L. * **Deflection under the load** at :math:`x = L/2`: .. math:: :label: propped-cantilever-pt-deflection v\bigl(\tfrac{L}{2}\bigr) \;=\; -\,\frac{7\, P\, L^{3}}{768\, E\, I}. The first three are extracted from the BEAM2 reaction-force output; the fourth comes from the nodal displacement at the load point. Comparison to the cantilever with the same midspan load (no prop): removing the prop reaction takes the deflection at :math:`x = L/2` from :math:`-7 P L^{3} / 768 EI` to :math:`-P L^{3}/(24 EI) = -32 P L^{3} / 768 EI` — a factor of about :math:`32/7 \approx 4.6\times` larger. The prop kills most of the bending compliance even though it sits an entire half-span away from the load. Implementation -------------- A 40-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam. The root is fully clamped (every DOF pinned); the tip is "simply supported" (only :math:`u_y = 0`). Out-of-plane DOFs are pinned across the line so the response stays strictly 2D in the :math:`x`-:math:`y` plane. Hermite-cubic shape functions interpolate Bernoulli kinematics exactly on the polynomial moment field, so the FE response under this single concentrated load matches the closed form to machine precision at every node. References ---------- * Timoshenko, S. P. (1955) *Strength of Materials, Part I: Elementary Theory and Problems*, 3rd ed., Van Nostrand, §17 (force-method solution of 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 13a (propped cantilever, concentrated load at any point — midspan special case). * Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley, §16.2 (statically-indeterminate beam FE). Vendor cross-references ----------------------- ====================================== ========================= ============================================ Source Reported δ_mid [m] Problem ID / location ====================================== ========================= ============================================ Closed form (Gere & Goodno) 8.75 × 10⁻⁵ §10.3 Table 10-1 Case 6 Timoshenko (1955) 8.75 × 10⁻⁵ SoM Part I §5.8 MAPDL Verification Manual 8.75 × 10⁻⁵ VM-41 family (propped cantilever) Abaqus Verification Manual 8.75 × 10⁻⁵ AVM 1.5.x propped-cantilever family ====================================== ========================= ============================================ """ from __future__ import annotations import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS # %% # Problem data # ------------ 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] P = 1.0e3 # transverse load at midspan [N, downward = -y] EI = E * I_z # Closed-form references — Roark Table 8 case 13a (midspan special case). R_root_pub = 11.0 * P / 16.0 R_prop_pub = 5.0 * P / 16.0 M_root_pub = -3.0 * P * L / 16.0 v_midspan_pub = -7.0 * P * L**3 / (768.0 * EI) # Sanity-check the unpropped cantilever for the comparison print. v_midspan_unpropped = -P * L**3 / (24.0 * EI) print("Propped cantilever — central point load") print(f" L = {L} m, P = {P} N at x = L/2 = {L / 2} m") print(f" E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2") print() print("Closed-form references (Roark Table 8 case 13a):") print(f" R_prop = {R_prop_pub:+.4e} N (= 5 P / 16)") print(f" R_root = {R_root_pub:+.4e} N (= 11 P / 16)") print(f" M_root = {M_root_pub:+.4e} N m (= -3 P L / 16)") print(f" v(L/2) = {v_midspan_pub * 1e3:+.4e} mm (= -7 P L^3 / 768 EI)") print() print(f" v(L/2) unpropped cantilever for comparison = {v_midspan_unpropped * 1e3:+.4e} mm") print( f" → propping the tip stiffens midspan deflection by " f"{abs(v_midspan_unpropped / v_midspan_pub):.2f}×" ) # %% # 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 # ------------------- # # Root at x=0: full clamp (UX=UY=UZ=ROTX=ROTY=ROTZ = 0). # Prop at x=L: simply supported (only UY=0; in-plane rotation free). # Out-of-plane motion (UZ, ROTX, ROTY) suppressed across the line so # the response stays strictly 2D. m.fix(nodes=1, dof="ALL") m.fix(nodes=int(N_ELEM + 1), dof="UY") for i in range(N_ELEM + 1): m.fix(nodes=int(i + 1), dof="UZ") m.fix(nodes=int(i + 1), dof="ROTX") m.fix(nodes=int(i + 1), dof="ROTY") # %% # Load at midspan # --------------- i_mid = N_ELEM // 2 m.apply_force(int(i_mid + 1), fy=-P) # %% # Solve + extract midspan deflection # ---------------------------------- 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_midspan_fe = _node_dof(int(i_mid + 1), 1) err_midspan = (abs(v_midspan_fe) - abs(v_midspan_pub)) / abs(v_midspan_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(L/2) midspan':<22} {v_midspan_fe * 1e3:>10.4f} mm " f"{v_midspan_pub * 1e3:>10.4f} mm {err_midspan:>+8.4f}%" ) assert abs(err_midspan) < 0.1, f"v(L/2) deviation {err_midspan:.4f} % exceeds 0.1 %" # %% # Plot the deflection profile against the closed form # --------------------------------------------------- # # The closed-form deflection splits at the load point. For # :math:`0 \le x \le L/2`, # :math:`E I\, v(x) = -3 P L x^{2} / 32 + 11 P x^{3} / 96`. For # :math:`L/2 \le x \le L`, # :math:`E I\, v(x) = 5 P L x^{2} / 64 - (5 P / 96)(x - L/2)^{3} # - 11 P L^{2} x / 128 + 11 P L^{3} / 768`. def v_closed_form(x: float) -> float: if x <= L / 2: return (-3.0 * P * L * x**2 / 32.0 + 11.0 * P * x**3 / 96.0) / EI * (-1.0) return ( ( 5.0 * P * L * x**2 / 64.0 - (5.0 * P / 96.0) * (x - L / 2) ** 3 - 11.0 * P * L**2 * x / 128.0 + 11.0 * P * L**3 / 768.0 ) / EI * (-1.0) ) fe_uy = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)]) 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 13a)", ) ax.plot(xs, fe_uy * 1e3, "o", color="#d62728", markersize=5, label="BEAM2 FE") ax.axvline(L / 2, color="grey", ls="--", lw=1.0, label=f"load at x = L/2 = {L / 2:.3f} m") ax.axhline(0.0, color="black", lw=0.5) ax.set_xlabel("x [m]") ax.set_ylabel("v(x) [mm] (downward deflection negative)") ax.set_title("Propped cantilever, midspan load — FE vs closed form", fontsize=11) ax.legend(loc="lower right", fontsize=9) ax.grid(True, ls=":", alpha=0.5) fig.tight_layout() fig.show() # %% # Render the deformed beam # ------------------------ 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_midspan_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([[L, 0.0, 0.0]]), render_points_as_spheres=True, point_size=18, color="#888888", label="prop (UY = 0)", ) plotter.add_points( np.array([[L / 2, scale * v_midspan_fe, 0.0]]), render_points_as_spheres=True, point_size=14, color="#d62728", label=f"load — v(L/2) = {v_midspan_fe * 1e3:.4f} mm", ) plotter.view_xy() plotter.add_legend() plotter.show() # %% # Take-aways # ---------- print() print("Take-aways:") print(f" • v(L/2) within {abs(err_midspan):.4f} % of -7 P L³ / (768 E I).") print( f" • Adding the prop at x = L stiffens the beam against this midspan load by " f"{abs(v_midspan_unpropped / v_midspan_pub):.2f}× compared to a free-tip cantilever." )