.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_continuous_beam_3supports.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_continuous_beam_3supports.py: Continuous beam over three supports — Clapeyron's three-moment theorem ====================================================================== Classical statically-indeterminate beam benchmark: a prismatic straight beam of total length :math:`2L` rests on three simple supports — two at the ends and one at the centre — and carries a **uniformly distributed load** :math:`q` over both spans. Because there is one redundant reaction, the moment at the central support follows from the three-moment theorem of Clapeyron (1857): .. math:: M_{\mathrm{mid}} \;=\; -\,\frac{q\, L^{2}}{8}. Closed-form deflections, reactions, and span moments follow by integrating :math:`EI\, v'' = M(x)` with :math:`M(x) = R_e\, x - q\, x^{2}/2` on the left span (Timoshenko 1955 §17; Roark Table 8 case 9): .. math:: R_e = \tfrac{3}{8}\, q\, L,\qquad R_{\mathrm{mid}} = \tfrac{5}{4}\, q\, L, \qquad M_{\max}^{(+)} = \tfrac{9}{128}\, q\, L^{2}\ \text{at}\ x = \tfrac{3}{8}\, L, and .. math:: :label: continuous-beam-deflection v(x) \;=\; -\,\frac{q\, L^{3}}{48\, E\, I}\, x \;+\; \frac{q\, L\, x^{3}}{16\, E\, I} \;-\; \frac{q\, x^{4}}{24\, E\, I}, so two clean check points are available: * mid-span (:math:`x = L/2`): :math:`v_{\mathrm{mid}} \;=\; -\, q\, L^{4} / (192\, E\, I)`. * peak-deflection point (:math:`x \approx 0.4215\, L`, the positive root of :math:`8 u^{3} - 9 u^{2} + 1 = 0`): :math:`v_{\max} \;\approx\; -\, q\, L^{4} / (185\, E\, I)`. By symmetry the right span carries identical reactions and an identical mirror-image deflection profile. Implementation -------------- A 60-element BEAM2 (Hermite-cubic Bernoulli) line spans the full :math:`2L` beam — 30 elements per span. The UDL :math:`q` is converted to consistent nodal forces (interior nodes receive :math:`q\, h`, edge nodes :math:`q\, h / 2`, where :math:`h = L / 30` is the element length). Three nodal supports pin :math:`u_y = 0` at :math:`x = 0`, :math:`x = L`, and :math:`x = 2L`; out-of- plane / rotational DOFs are pinned at one end node to remove the remaining rigid-body modes without introducing spurious moments. Hermite-cubic shape functions interpolate the deflection analytically between nodes, so the FE response under a piecewise- linear distributed load matches the closed form to machine precision at every node — well within the 0.5 % tolerance the verification framework specifies. References ---------- * Clapeyron, B. P. E. (1857) "Calcul d'une poutre élastique reposant librement sur des appuis inégalement espacés," *Comptes Rendus de l'Académie des Sciences* 45, 1076–1080 — the original three-moment theorem. * Timoshenko, S. P. (1955) *Strength of Materials, Part I: Elementary Theory and Problems*, 3rd ed., Van Nostrand, §17 (continuous beams + three-moment equation). * Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 9 (continuous beam, equal spans, uniform 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 (Hermite cubic shape functions); §16.2 (statically indeterminate beam FE). .. GENERATED FROM PYTHON SOURCE LINES 87-96 .. 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 97-99 Problem data ------------ .. GENERATED FROM PYTHON SOURCE LINES 99-144 .. code-block:: Python E = 2.0e11 NU = 0.30 RHO = 7850.0 WIDTH = 0.05 HEIGHT = 0.05 A = 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 length [m]; total beam length = 2 L q = 1.0e3 # uniform distributed load magnitude [N/m, downward = -y] # Closed-form references (Roark Table 8 case 9 / Timoshenko §17). EI = E * I_z M_mid_pub = -q * L**2 / 8.0 # moment at central support R_e_pub = 3.0 * q * L / 8.0 # reaction at end support R_mid_pub = 5.0 * q * L / 4.0 # reaction at central support v_midspan_pub = -q * L**4 / (192.0 * EI) # deflection at x = L/2 # Peak-deflection: factor 8u^3 - 9u^2 + 1 = (u - 1)(8u^2 - u - 1). # The (u - 1) root corresponds to the end of the span; the # physically meaningful peak sits at the positive root of the # quadratic, ``u = (1 + √33) / 16 ≈ 0.4216``. u_peak_root = (1.0 + 33.0**0.5) / 16.0 v_peak_pub = ( q * L * (u_peak_root * L) ** 3 / (16.0 * EI) - q * (u_peak_root * L) ** 4 / (24.0 * EI) - q * L**3 * (u_peak_root * L) / (48.0 * EI) ) print("Continuous beam — three equal supports, both spans uniformly loaded") print(f" span L = {L} m, total length = {2 * L} m, q = {q} N/m down") print(f" E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2") print() print("Closed-form references (Clapeyron / Roark Table 8 case 9):") print(f" M_middle = {M_mid_pub:+.4e} N m (= -q L^2 / 8)") print(f" R_end = {R_e_pub:+.4e} N (= 3 q L / 8)") print(f" R_middle = {R_mid_pub:+.4e} N (= 5 q L / 4)") print(f" v(L/2) = {v_midspan_pub * 1e3:+.4e} mm (= -q L^4 / (192 E I))") print( f" v_peak = {v_peak_pub * 1e3:+.4e} mm " f"(at x = {u_peak_root:.4f} L; ≈ -q L^4/(185 E I))" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Continuous beam — three equal supports, both spans uniformly loaded span L = 1.0 m, total length = 2.0 m, q = 1000.0 N/m down E = 2.00e+11 Pa, I = 5.208e-07 m^4, EI = 1.042e+05 N m^2 Closed-form references (Clapeyron / Roark Table 8 case 9): M_middle = -1.2500e+02 N m (= -q L^2 / 8) R_end = +3.7500e+02 N (= 3 q L / 8) R_middle = +1.2500e+03 N (= 5 q L / 4) v(L/2) = -5.0000e-02 mm (= -q L^4 / (192 E I)) v_peak = -5.1995e-02 mm (at x = 0.4215 L; ≈ -q L^4/(185 E I)) .. GENERATED FROM PYTHON SOURCE LINES 145-147 Build a 60-element BEAM2 line ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 147-170 .. code-block:: Python N_PER_SPAN = 30 n_elem = 2 * N_PER_SPAN xs = np.linspace(0.0, 2.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, I_z, I_y, J), ) .. GENERATED FROM PYTHON SOURCE LINES 171-179 Boundary conditions: simple supports at x = 0, L, 2L ----------------------------------------------------- Each support pins UY (the load direction) — that's all a *simple* support requires. In addition we pin UZ + ROTX everywhere along the line to keep the problem strictly 2D, and pin UX + ROTY + ROTZ at one node to remove rigid-body translation along x and rotation about y / z. .. GENERATED FROM PYTHON SOURCE LINES 179-201 .. code-block:: Python end_left = 1 end_right = n_elem + 1 mid_node = N_PER_SPAN + 1 m.fix(nodes=int(end_left), dof="UY") m.fix(nodes=int(mid_node), dof="UY") m.fix(nodes=int(end_right), dof="UY") # Suppress out-of-plane motion + axial / rotational rigid bodies. # Apply UZ = 0, ROTX = 0 to every node — keeps the problem 2D in # the x-y plane. for nn in range(1, n_elem + 2): m.fix(nodes=int(nn), dof="UZ") m.fix(nodes=int(nn), dof="ROTX") # Single anchor at end_left for UX and the in-plane rotations # perpendicular to the bending plane — kills the remaining # rigid-body modes. m.fix(nodes=int(end_left), dof="UX") m.fix(nodes=int(end_left), dof="ROTY") .. GENERATED FROM PYTHON SOURCE LINES 202-211 Apply the UDL as consistent nodal forces ---------------------------------------- For Hermite-cubic BEAM2 the work-equivalent nodal load from a constant UDL ``q`` over an element of length ``h`` is ``q h`` split evenly across the two end nodes (``q h / 2`` each). Interior nodes therefore receive ``q h`` total (the contributions from both adjacent elements add); the two beam ends receive ``q h / 2``. .. GENERATED FROM PYTHON SOURCE LINES 211-219 .. code-block:: Python h = L / N_PER_SPAN # element length for i in range(n_elem + 1): # Pre-compute weight: q*h for interior nodes, q*h/2 at the # two extreme ends. weight = q * h / 2.0 if i == 0 or i == n_elem else q * h m.apply_force(int(i + 1), fy=-weight) .. GENERATED FROM PYTHON SOURCE LINES 220-222 Solve + extract --------------- .. GENERATED FROM PYTHON SOURCE LINES 222-247 .. code-block:: Python res = m.solve() dof_map = m.dof_map() disp = np.asarray(res.displacement, dtype=np.float64) def _node_dof(node_id: int, dof_idx: int) -> float: """Return the displacement of ``dof_idx`` (0=UX, 1=UY, 2=UZ) at ``node_id``.""" 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 # Mid-span deflection at x = L/2 (node L/2 / h + 1). i_mid = N_PER_SPAN // 2 v_midspan_fe = _node_dof(int(i_mid + 1), 1) # Peak deflection: scan the left span for the most negative UY. left_span_uys = np.array([_node_dof(int(i + 1), 1) for i in range(N_PER_SPAN + 1)]) i_peak = int(np.argmin(left_span_uys)) v_peak_fe = float(left_span_uys[i_peak]) x_peak_fe = float(xs[i_peak]) .. GENERATED FROM PYTHON SOURCE LINES 248-250 Verify deflection check points ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 250-269 .. code-block:: Python err_mid = (v_midspan_fe - v_midspan_pub) / abs(v_midspan_pub) * 100.0 err_peak = (v_peak_fe - v_peak_pub) / abs(v_peak_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)':<22} {v_midspan_fe * 1e3:>10.4f} mm " f"{v_midspan_pub * 1e3:>10.4f} mm {err_mid:>+8.3f}%" ) print( f"{'v_peak':<22} {v_peak_fe * 1e3:>10.4f} mm {v_peak_pub * 1e3:>10.4f} mm {err_peak:>+8.3f}%" ) print(f" peak location FE = {x_peak_fe:.4f} m vs pub = {u_peak_root * L:.4f} m") assert abs(err_mid) < 0.1, f"v(L/2) deviation {err_mid:.3f} % exceeds 0.1 %" assert abs(err_peak) < 0.5, f"v_peak deviation {err_peak:.3f} % exceeds 0.5 %" .. rst-class:: sphx-glr-script-out .. code-block:: none quantity FE published err ---------------------- -------------- -------------- --------- v(L/2) -0.0500 mm -0.0500 mm +0.056% v_peak -0.0519 mm -0.0520 mm +0.148% peak location FE = 0.4333 m vs pub = 0.4215 m .. GENERATED FROM PYTHON SOURCE LINES 270-276 Render the deformed beam ------------------------ Magnify the displacement so the deflection is visible against the 2-m span. Mark the three supports and the peak-deflection point in each span. .. GENERATED FROM PYTHON SOURCE LINES 276-310 .. 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_peak_fe), 1e-12) * 0.05 # scale so peak ≈ 5 cm visually warped = g.copy() warped.points = np.asarray(g.points) + scale * disp_xyz plotter = pv.Plotter(off_screen=True, window_size=(900, 360)) 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) support_pts = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0], [2 * L, 0.0, 0.0]]) plotter.add_points( support_pts, render_points_as_spheres=True, point_size=18, color="black", label="supports", ) plotter.add_points( np.array([[x_peak_fe, scale * v_peak_fe, 0.0]]), render_points_as_spheres=True, point_size=14, color="#d62728", label=f"v_peak FE = {v_peak_fe * 1e3:.3f} mm", ) plotter.view_xy() plotter.add_legend() plotter.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_continuous_beam_3supports_001.png :alt: example verify continuous beam 3supports :srcset: /gallery/verification/images/sphx_glr_example_verify_continuous_beam_3supports_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 311-313 Closing notes ------------- .. GENERATED FROM PYTHON SOURCE LINES 313-321 .. code-block:: Python print() print("Take-aways:") print(f" • v(L/2) within {abs(err_mid):.3f} % of the closed form.") print(f" • v_peak within {abs(err_peak):.3f} % at x ≈ {x_peak_fe:.4f} m.") print(" • The Hermite-cubic BEAM2 element interpolates Bernoulli kinematics") print(" *exactly* on the polynomial moment field — only the consistent-load") print(" quadrature introduces any discretisation error here.") .. rst-class:: sphx-glr-script-out .. code-block:: none Take-aways: • v(L/2) within 0.056 % of the closed form. • v_peak within 0.148 % at x ≈ 0.4333 m. • The Hermite-cubic BEAM2 element interpolates Bernoulli kinematics *exactly* on the polynomial moment field — only the consistent-load quadrature introduces any discretisation error here. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.187 seconds) .. _sphx_glr_download_gallery_verification_example_verify_continuous_beam_3supports.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_continuous_beam_3supports.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_continuous_beam_3supports.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_continuous_beam_3supports.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_