.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_cantilever_tip_moment.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_tip_moment.py: Cantilever beam under a tip moment ================================== Pure-bending end-loaded prismatic beam, the simplest test of moment-driven Bernoulli kinematics. A clamped-free cantilever of length :math:`L` carries a moment :math:`M_0` about the out-of-plane axis applied at its tip; equilibrium and elementary integration give a closed form (Roark & Young 1989, Table 8 case 4; Timoshenko 1955 §32): .. math:: \theta(x) = \frac{M_0\, x}{E I}, \qquad v(x) = \frac{M_0\, x^{2}}{2\, E I}, so the slope and tip rotation / deflection are .. math:: \theta_\mathrm{tip} = \frac{M_0\, L}{E I}, \qquad v_\mathrm{tip} = \frac{M_0\, L^{2}}{2\, E I}, and the curvature :math:`\kappa = M_0 / (E I)` is **constant** along the entire span — a cantilever under pure end moment deforms into a circular arc of radius :math:`E I / M_0`. Implementation -------------- A 10-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam. Boundary conditions: clamp all 6 DOFs at the root. Loading: a single :math:`+M_z` moment at the tip node. Hermite cubic shape functions interpolate displacement and slope simultaneously, so the FE response under a pure tip moment matches the analytic parabolic deflection profile to machine precision at every node. References ---------- * Timoshenko, S. P. (1955) *Strength of Materials, Part I*, 3rd ed., Van Nostrand, §32. * Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 4 (cantilever, end moment). * 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 52-62 .. code-block:: Python from __future__ import annotations import numpy as np import pyvista as pv from vtkmodules.util.vtkConstants import VTK_LINE import femorph_solver from femorph_solver import ELEMENTS .. GENERATED FROM PYTHON SOURCE LINES 63-65 Problem data ------------ .. GENERATED FROM PYTHON SOURCE LINES 65-91 .. code-block:: Python E = 2.0e11 # Pa (steel) NU = 0.30 RHO = 7850.0 # kg/m^3 b = 0.050 # square cross-section side [m] A = b * b I = b**4 / 12.0 # noqa: E741 J = 2.0 * I L = 1.0 # span [m] M0 = 1.0e3 # tip moment [N m] about +z N_ELEM = 10 # 10 segments → exact at every node for Bernoulli kinematics # Closed-form quantities ----------------------------------------------- theta_tip_published = M0 * L / (E * I) v_tip_published = M0 * L**2 / (2.0 * E * I) kappa_published = M0 / (E * I) R_arc = 1.0 / kappa_published # radius of the deformed arc print(f"M0 = {M0:.1f} N m, L = {L:.2f} m, EI = {E * I:.3e} N m^2") print(f"θ_tip = M L / (E I) = {theta_tip_published:.6e} rad") print(f"v_tip = M L^2 / (2 E I) = {v_tip_published:.6e} m") print(f"κ = M / (E I) = {kappa_published:.6e} 1/m (constant)") print(f"R_arc = E I / M = {R_arc:.4f} m (radius of deformed arc)") .. rst-class:: sphx-glr-script-out .. code-block:: none M0 = 1000.0 N m, L = 1.00 m, EI = 1.042e+05 N m^2 θ_tip = M L / (E I) = 9.600000e-03 rad v_tip = M L^2 / (2 E I) = 4.800000e-03 m κ = M / (E I) = 9.600000e-03 1/m (constant) R_arc = E I / M = 104.1667 m (radius of deformed arc) .. GENERATED FROM PYTHON SOURCE LINES 92-94 Build a 10-element BEAM2 line clamped at x=0 -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 94-116 .. code-block:: Python points = np.array( [[i * L / N_ELEM, 0.0, 0.0] for i in range(N_ELEM + 1)], dtype=np.float64, ) cells_list: list[int] = [] for i in range(N_ELEM): cells_list.extend([2, i, i + 1]) cells = np.asarray(cells_list, dtype=np.int64) cell_types = np.full(N_ELEM, VTK_LINE, dtype=np.uint8) grid = pv.UnstructuredGrid(cells, cell_types, points) m = femorph_solver.Model.from_grid(grid) m.assign( ELEMENTS.BEAM2, material={"EX": E, "PRXY": NU, "DENS": RHO}, real=(A, I, I, J), ) m.fix(nodes=[1], dof="ALL") # full clamp at root m.apply_force(N_ELEM + 1, mz=M0) # tip moment about +z .. GENERATED FROM PYTHON SOURCE LINES 117-119 Static solve + extract tip rotation / deflection ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 119-138 .. code-block:: Python res = m.solve() dof = m.dof_map() tip_uy_idx = np.where((dof[:, 0] == N_ELEM + 1) & (dof[:, 1] == 1))[0][0] tip_rotz_idx = np.where((dof[:, 0] == N_ELEM + 1) & (dof[:, 1] == 5))[0][0] v_tip = float(res.displacement[tip_uy_idx]) theta_tip = float(res.displacement[tip_rotz_idx]) err_v = (v_tip - v_tip_published) / v_tip_published err_th = (theta_tip - theta_tip_published) / theta_tip_published print() print("tip response (femorph-solver) → (analytical)") print(f" v_tip = {v_tip:+.6e} → {v_tip_published:+.6e} ({err_v * 100:+.6f}%)") print(f" θ_tip = {theta_tip:+.6e} → {theta_tip_published:+.6e} ({err_th * 100:+.6f}%)") assert abs(err_v) < 1e-8, f"tip deflection error {err_v:.2e} too large for Bernoulli BEAM2" assert abs(err_th) < 1e-8, f"tip rotation error {err_th:.2e} too large for Bernoulli BEAM2" .. rst-class:: sphx-glr-script-out .. code-block:: none tip response (femorph-solver) → (analytical) v_tip = +4.800000e-03 → +4.800000e-03 (-0.000000%) θ_tip = +9.600000e-03 → +9.600000e-03 (-0.000000%) .. GENERATED FROM PYTHON SOURCE LINES 139-144 Verify the parabolic deflection profile --------------------------------------- Every node should satisfy :math:`v(x) = M_0 x^2 / (2 E I)` to machine precision. .. GENERATED FROM PYTHON SOURCE LINES 144-155 .. code-block:: Python x_nodes = points[:, 0] v_published_profile = M0 * x_nodes**2 / (2.0 * E * I) v_computed_profile = np.zeros(N_ELEM + 1) for i in range(N_ELEM + 1): rows = np.where((dof[:, 0] == i + 1) & (dof[:, 1] == 1))[0] v_computed_profile[i] = float(res.displacement[rows[0]]) np.testing.assert_allclose(v_computed_profile, v_published_profile, atol=1e-14, rtol=1e-10) print() print("OK — full deflection profile v(x) = M_0 x^2 / (2 E I) matches at every node.") .. rst-class:: sphx-glr-script-out .. code-block:: none OK — full deflection profile v(x) = M_0 x^2 / (2 E I) matches at every node. .. GENERATED FROM PYTHON SOURCE LINES 156-163 Verify constant curvature κ = M / (EI) -------------------------------------- A second central difference on the slope :math:`\theta_z(x)` at every interior node should recover the constant :math:`\kappa = M_0 / (E I)` to within the Bernoulli-FE-at-node precision floor. .. GENERATED FROM PYTHON SOURCE LINES 163-176 .. code-block:: Python theta_profile = np.zeros(N_ELEM + 1) for i in range(N_ELEM + 1): rows = np.where((dof[:, 0] == i + 1) & (dof[:, 1] == 5))[0] theta_profile[i] = float(res.displacement[rows[0]]) kappa_estimates = np.diff(theta_profile) / np.diff(x_nodes) print( f"κ estimates (per element): mean = {kappa_estimates.mean():.6e}, " f"std = {kappa_estimates.std():.2e}" ) np.testing.assert_allclose(kappa_estimates, kappa_published, rtol=1e-10) print(f"OK — κ uniformly equals M_0/(E I) = {kappa_published:.6e} 1/m within 1e-10.") .. rst-class:: sphx-glr-script-out .. code-block:: none κ estimates (per element): mean = 9.600000e-03, std = 6.48e-16 OK — κ uniformly equals M_0/(E I) = 9.600000e-03 1/m within 1e-10. .. GENERATED FROM PYTHON SOURCE LINES 177-179 Render the deflected line + reference circular arc -------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 179-205 .. code-block:: Python x_dense = np.linspace(0, L, 200) v_dense = M0 * x_dense**2 / (2.0 * E * I) arc_pts = np.column_stack([x_dense, v_dense, np.zeros_like(x_dense)]) arc = pv.PolyData(arc_pts) warp = grid.copy() warp.points = grid.points + np.column_stack( [np.zeros(N_ELEM + 1), 200.0 * v_computed_profile, np.zeros(N_ELEM + 1)] ) warp["uy"] = v_computed_profile plotter = pv.Plotter(off_screen=True, window_size=(720, 320)) plotter.add_mesh(grid, color="grey", opacity=0.3, line_width=2) plotter.add_mesh( arc, style="points", color="black", opacity=0.4, point_size=3, label="closed form v(x)", ) plotter.add_mesh(warp, scalars="uy", line_width=5, cmap="viridis") plotter.view_xy() plotter.camera.zoom(1.05) plotter.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_cantilever_tip_moment_001.png :alt: example verify cantilever tip moment :srcset: /gallery/verification/images/sphx_glr_example_verify_cantilever_tip_moment_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.159 seconds) .. _sphx_glr_download_gallery_verification_example_verify_cantilever_tip_moment.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_tip_moment.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_cantilever_tip_moment.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_cantilever_tip_moment.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_