.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_cantilever_midspan_load.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_verification_example_verify_cantilever_midspan_load.py: 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). .. GENERATED FROM PYTHON SOURCE LINES 72-81 .. code-block:: Python from __future__ import annotations import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS .. GENERATED FROM PYTHON SOURCE LINES 82-84 Problem data — slender beam, midspan load ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 84-129 .. code-block:: Python 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") .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 130-132 Build a 40-element BEAM2 line ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 132-154 .. code-block:: Python 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), ) .. GENERATED FROM PYTHON SOURCE LINES 155-157 Boundary conditions: full clamp at x = 0 ---------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 157-160 .. code-block:: Python m.fix(nodes=1, dof="ALL") .. GENERATED FROM PYTHON SOURCE LINES 161-163 Apply the load at the closest node to x = a ------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 163-167 .. code-block:: Python i_load = int(np.argmin(np.abs(xs - a))) m.apply_force(int(i_load + 1), fy=-P) .. GENERATED FROM PYTHON SOURCE LINES 168-170 Suppress out-of-plane motion at every node so the strictly-2D response matches the closed-form Bernoulli kinematics. .. GENERATED FROM PYTHON SOURCE LINES 170-176 .. code-block:: Python 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") .. GENERATED FROM PYTHON SOURCE LINES 177-179 Solve + extract --------------- .. GENERATED FROM PYTHON SOURCE LINES 179-222 .. code-block:: Python 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 %" .. rst-class:: sphx-glr-script-out .. code-block:: none 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% .. GENERATED FROM PYTHON SOURCE LINES 223-225 Plot the deflection profile against the closed form --------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 225-251 .. code-block:: Python 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() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_cantilever_midspan_load_001.png :alt: Cantilever with off-tip point load — FE vs closed form :srcset: /gallery/verification/images/sphx_glr_example_verify_cantilever_midspan_load_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 252-254 Render the deformed shape ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 254-287 .. code-block:: Python 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() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_cantilever_midspan_load_002.png :alt: example verify cantilever midspan load :srcset: /gallery/verification/images/sphx_glr_example_verify_cantilever_midspan_load_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 288-290 Take-aways ---------- .. GENERATED FROM PYTHON SOURCE LINES 290-301 .. code-block:: Python 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)." ) .. rst-class:: sphx-glr-script-out .. code-block:: none 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). .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.283 seconds) .. _sphx_glr_download_gallery_verification_example_verify_cantilever_midspan_load.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_verify_cantilever_midspan_load.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_cantilever_midspan_load.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_cantilever_midspan_load.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_