.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_ss_beam_central_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_ss_beam_central_load.py: Simply-supported beam under a central point load ================================================ Classical statically-determinate first-order beam: simple supports (pin + roller) at :math:`x = 0` and :math:`x = L`, with a single transverse point load :math:`P` at midspan. Symmetry plus equilibrium gives the elementary closed form .. math:: R_\mathrm{left} = R_\mathrm{right} = \tfrac{P}{2}, \qquad M_\mathrm{max} = \tfrac{P L}{4}\;\text{at}\;x = L/2, with mid-span and quarter-span deflections (Roark & Young 1989, Table 8 case 1; Timoshenko 1955 §32) .. math:: \delta(L/2) = \frac{P L^{3}}{48\, E I}, \qquad \delta(L/4) = \frac{11\, P L^{3}}{768\, E I}. Implementation -------------- A 20-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam. Boundary conditions: * node 1 (pin): UY = UZ = UX = 0; ROTX, ROTY suppressed; ROTZ free (slope rotates). * node :math:`N_e + 1` (roller): UY = UZ = 0; ROTX, ROTY suppressed; UX and ROTZ free. * tip load :math:`-P\,\hat y` at the central node :math:`N_e/2 + 1`. The Hermite-cubic basis recovers Euler-Bernoulli kinematics exactly under nodal point loads, so the FE solution at every node equals the analytical value to machine precision. References ---------- * Timoshenko, S. P. (1955) *Strength of Materials, Part I*, 3rd ed., Van Nostrand, §32 (concentrated transverse loads). * Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 1 (simply-supported beam, central point load). * 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. .. GENERATED FROM PYTHON SOURCE LINES 53-63 .. 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 64-66 Problem data ------------ .. GENERATED FROM PYTHON SOURCE LINES 66-93 .. 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] P = 5.0e3 # tip load magnitude [N], applied downward (-y) N_ELEM = 20 # must be even so a node lands at midspan # Closed-form quantities ----------------------------------------------- R_left_published = P / 2.0 R_right_published = P / 2.0 delta_mid_published = P * L**3 / (48.0 * E * I) delta_quarter_published = 11.0 * P * L**3 / (768.0 * E * I) M_max_published = P * L / 4.0 print(f"P = {P:.1f} N, L = {L:.2f} m, EI = {E * I:.3e} N m^2") print(f"R_L = R_R = P/2 = {R_left_published:.4f} N") print(f"M_max = P L / 4 = {M_max_published:.4f} N m at midspan") print(f"δ(L/2) = P L^3/(48 EI) = {delta_mid_published:.4e} m") print(f"δ(L/4) = 11 P L^3/(768EI) = {delta_quarter_published:.4e} m") .. rst-class:: sphx-glr-script-out .. code-block:: none P = 5000.0 N, L = 1.00 m, EI = 1.042e+05 N m^2 R_L = R_R = P/2 = 2500.0000 N M_max = P L / 4 = 1250.0000 N m at midspan δ(L/2) = P L^3/(48 EI) = 1.0000e-03 m δ(L/4) = 11 P L^3/(768EI) = 6.8750e-04 m .. GENERATED FROM PYTHON SOURCE LINES 94-96 Build a 20-element BEAM2 line ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 96-115 .. 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), ) .. GENERATED FROM PYTHON SOURCE LINES 116-123 Boundary conditions ------------------- Standard simply-supported beam: * left pin — UY, UZ, UX = 0; ROTX, ROTY suppressed; ROTZ free. * right roller — UY, UZ = 0; ROTX, ROTY suppressed; UX & ROTZ free. .. GENERATED FROM PYTHON SOURCE LINES 123-138 .. code-block:: Python # Left pin — pin axial too so the system is statically determinate. m.fix(nodes=[1], dof="UX") m.fix(nodes=[1], dof="UY") m.fix(nodes=[1], dof="UZ") m.fix(nodes=[1], dof="ROTX") m.fix(nodes=[1], dof="ROTY") # Right roller — only UY pinned (vertical reaction); UX free. right = N_ELEM + 1 m.fix(nodes=[right], dof="UY") m.fix(nodes=[right], dof="UZ") m.fix(nodes=[right], dof="ROTX") m.fix(nodes=[right], dof="ROTY") .. GENERATED FROM PYTHON SOURCE LINES 139-146 Apply the central point load ---------------------------- A single :math:`-P\hat y` force at the midspan node. No distributed-load consistent vector needed — the Hermite cubic basis captures the response of a point load on a Bernoulli beam exactly at every node. .. GENERATED FROM PYTHON SOURCE LINES 146-150 .. code-block:: Python mid_node = N_ELEM // 2 + 1 m.apply_force(mid_node, fy=-P) .. GENERATED FROM PYTHON SOURCE LINES 151-153 Static solve + reaction extraction ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 153-176 .. code-block:: Python res = m.solve() dof = m.dof_map() def _react(node: int, dof_idx: int) -> float: rows = np.where((dof[:, 0] == node) & (dof[:, 1] == dof_idx))[0] return float(res.reaction[rows[0]]) if len(rows) else 0.0 R_left = _react(1, 1) R_right = _react(right, 1) print() print("reactions (femorph-solver) → (analytical)") print(f" R_left_UY = {R_left:+.4f} → {+R_left_published:+.4f} N") print(f" R_right_UY = {R_right:+.4f} → {+R_right_published:+.4f} N") # Vertical equilibrium and exact agreement. assert np.isclose(R_left + R_right, P, rtol=1e-12), "vertical equilibrium failed" assert np.isclose(R_left, R_left_published, rtol=1e-12), "left reaction off" assert np.isclose(R_right, R_right_published, rtol=1e-12), "right reaction off" .. rst-class:: sphx-glr-script-out .. code-block:: none reactions (femorph-solver) → (analytical) R_left_UY = +2500.0000 → +2500.0000 N R_right_UY = +2500.0000 → +2500.0000 N .. GENERATED FROM PYTHON SOURCE LINES 177-179 Mid-span and quarter-span deflection ------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 179-197 .. code-block:: Python mid_uy = float(res.displacement[np.where((dof[:, 0] == mid_node) & (dof[:, 1] == 1))[0][0]]) err_mid = (mid_uy - (-delta_mid_published)) / (-delta_mid_published) print() print(f"δ(L/2) computed = {mid_uy:+.4e} m") print(f"δ(L/2) published = {-delta_mid_published:+.4e} m") print(f"relative error = {err_mid * 100:+.6f} %") assert abs(err_mid) < 1e-8, "mid-span deflection error too large for Bernoulli BEAM2" quarter_node = N_ELEM // 4 + 1 quarter_uy = float(res.displacement[np.where((dof[:, 0] == quarter_node) & (dof[:, 1] == 1))[0][0]]) err_quarter = (quarter_uy - (-delta_quarter_published)) / (-delta_quarter_published) print() print(f"δ(L/4) computed = {quarter_uy:+.4e} m") print(f"δ(L/4) published = {-delta_quarter_published:+.4e} m") print(f"relative error = {err_quarter * 100:+.6f} %") assert abs(err_quarter) < 1e-8, "quarter-span deflection error too large for Bernoulli BEAM2" .. rst-class:: sphx-glr-script-out .. code-block:: none δ(L/2) computed = -1.0000e-03 m δ(L/2) published = -1.0000e-03 m relative error = -0.000000 % δ(L/4) computed = -6.8750e-04 m δ(L/4) published = -6.8750e-04 m relative error = -0.000000 % .. GENERATED FROM PYTHON SOURCE LINES 198-200 Render the deflected line ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 200-230 .. code-block:: Python pts = np.asarray(grid.points) disp_y = np.zeros(N_ELEM + 1) for i in range(N_ELEM + 1): rows = np.where((dof[:, 0] == i + 1) & (dof[:, 1] == 1))[0] if len(rows): disp_y[i] = float(res.displacement[rows[0]]) warped = grid.copy() warped.points = pts + np.column_stack([np.zeros(N_ELEM + 1), 200.0 * disp_y, np.zeros(N_ELEM + 1)]) warped["uy"] = disp_y plotter = pv.Plotter(off_screen=True, window_size=(720, 280)) plotter.add_mesh(grid, color="grey", line_width=2, opacity=0.5) plotter.add_mesh(warped, scalars="uy", line_width=5, cmap="viridis") plotter.add_points( np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]]), render_points_as_spheres=True, point_size=18, color="black", ) plotter.add_points( np.array([[0.5 * L, 0.0, 0.0]]), render_points_as_spheres=True, point_size=14, color="#d62728", label=f"P = {P:.0f} N", ) plotter.view_xy() plotter.camera.zoom(1.05) plotter.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_ss_beam_central_load_001.png :alt: example verify ss beam central load :srcset: /gallery/verification/images/sphx_glr_example_verify_ss_beam_central_load_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.175 seconds) .. _sphx_glr_download_gallery_verification_example_verify_ss_beam_central_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_ss_beam_central_load.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_ss_beam_central_load.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_ss_beam_central_load.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_