r""" Cantilever with off-tip point load — load-position shape function ================================================================= A clamped-free prismatic beam of length :math:`L` with a single transverse point load :math:`P` applied at an interior position :math:`x = a \le L`. The Euler-Bernoulli closed form (Timoshenko 1955 §32; Roark Table 8 case 1) splits at the load: .. math:: v(x) = \begin{cases} \dfrac{P\, x^{2}\,(3 a - x)}{6\, E\, I}, & 0 \le x \le a \\[8pt] \dfrac{P\, a^{2}\,(3 x - a)}{6\, E\, I}, & a \le x \le L \end{cases} so the deflection grows as :math:`x^{2}\,(3 a - x)` while the load is "above" :math:`x` and as :math:`a^{2}\,(3 x - a)` past the load (a straight ramp of slope :math:`P\,a^{2} / (2\, E\, I)` away from the load). Three closed-form check points fall out: * **Deflection at the load** .. math:: v(a) = \frac{P\, a^{3}}{3\, E\, I}. * **Tip deflection** :math:`v(L)` .. math:: v(L) = \frac{P\, a^{2}\,(3 L - a)}{6\, E\, I}. * **Tip slope** :math:`\theta(L)` (= the slope past the load) .. math:: \theta(L) = \frac{P\, a^{2}}{2\, E\, I}. For the tip-load limit :math:`a = L` the formula collapses to :math:`v(L) = P L^{3} / (3 E I)`, recovering the :ref:`sphx_glr_gallery_verification_example_verify_cantilever_eb.py` result. For :math:`a = L/2` (midspan loading) the tip deflection is :math:`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 :math:`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") # %% # 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 # ------------------------------------------- i_load = int(np.argmin(np.abs(xs - a))) m.apply_force(int(i_load + 1), fy=-P) # %% # Suppress out-of-plane motion at every node so the strictly-2D # response matches the closed-form Bernoulli kinematics. 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") # %% # 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 %" # %% # 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)." )