.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_cantilever_torsion.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_torsion.py: Cantilever Saint-Venant torsion — rectangular cross-section =========================================================== Classical Saint-Venant elasticity verification: a clamped-free straight beam carrying an end torque :math:`T` at its free tip develops a uniform twist rate. The total tip rotation is (Roark & Young 1989 Table 20 case 4; Timoshenko & Goodier 1970 §90): .. math:: \varphi_\mathrm{tip} = \frac{T\, L}{G\, K}, \qquad G = \frac{E}{2 (1 + \nu)}, with :math:`K` the Saint-Venant torsion constant. For a **rectangular** cross-section :math:`b \times t` (with :math:`b \ge t`) the converged Saint-Venant series gives .. math:: K = \beta\bigl(b / t\bigr)\, b\, t^{3}, \qquad \beta(1) = 0.1406, (Roark §10 Table 20). Note: the polar moment :math:`J = b t (b^{2} + t^{2}) / 12` is **not** the torsion constant for non- circular sections — using :math:`J` instead of :math:`K` overstates the stiffness of a square shaft by ~12 %. Implementation -------------- Drives the existing :class:`~femorph_solver.validation.problems.CantileverTorsion` problem class on a 20-element BEAM2 (Cook Table 16.3-1 Hermite- cubic + Saint-Venant torsion) line. The validation suite uses a 2 % engineering tolerance — the tip twist matches the closed form to machine precision. References ---------- * Saint-Venant, B. (1855) "Mémoire sur la torsion des prismes", *Mém. Acad. Sci.* 14, 233–560. * Timoshenko, S. P. and Goodier, J. N. (1970) *Theory of Elasticity*, 3rd ed., McGraw-Hill, §90 (rectangular Saint-Venant torsion). * Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for Stress and Strain*, 6th ed., McGraw-Hill, Table 20 case 4 (rectangular Saint-Venant torsion constant). * Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley, Table 16.3-1 (BEAM2 Saint-Venant kernel). .. GENERATED FROM PYTHON SOURCE LINES 57-65 .. code-block:: Python from __future__ import annotations import numpy as np import pyvista as pv from femorph_solver.validation.problems import CantileverTorsion .. GENERATED FROM PYTHON SOURCE LINES 66-68 Build the model from the validation problem class -------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 68-72 .. code-block:: Python problem = CantileverTorsion() m = problem.build_model() .. GENERATED FROM PYTHON SOURCE LINES 73-75 Closed-form quantities ---------------------- .. GENERATED FROM PYTHON SOURCE LINES 75-87 .. code-block:: Python beta = 0.1406 # square cross-section (b/t = 1) G = problem.E / (2.0 * (1.0 + problem.nu)) K = beta * problem.width * problem.height**3 phi_published = problem.T * problem.L / (G * K) print(f"L = {problem.L} m, b × t = {problem.width} × {problem.height} m, T = {problem.T} N m") print(f"E = {problem.E / 1e9:.0f} GPa, ν = {problem.nu}, G = {G / 1e9:.2f} GPa") print(f"β(b/t = 1) = {beta} (Roark Table 20 case 4)") print(f"K = β·b·t^3 = {K:.4e} m^4 (Saint-Venant torsion constant)") print(f"φ_tip = T L / (G K) = {phi_published:.6e} rad ({np.degrees(phi_published):.4f} °)") .. rst-class:: sphx-glr-script-out .. code-block:: none L = 1.0 m, b × t = 0.05 × 0.05 m, T = 10.0 N m E = 200 GPa, ν = 0.3, G = 76.92 GPa β(b/t = 1) = 0.1406 (Roark Table 20 case 4) K = β·b·t^3 = 8.7875e-07 m^4 (Saint-Venant torsion constant) φ_tip = T L / (G K) = 1.479374e-04 rad (0.0085 °) .. GENERATED FROM PYTHON SOURCE LINES 88-90 Static solve + tip-twist extraction ----------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 90-104 .. code-block:: Python res = m.solve() phi_computed = problem.extract(m, res, "tip_twist") err_pct = (phi_computed - phi_published) / phi_published * 100.0 print() print(f"φ_tip computed = {phi_computed:.6e} rad") print(f"φ_tip published = {phi_published:.6e} rad") print(f"relative error = {err_pct:+.4f} %") # 2% tolerance — the validation class's tight engineering bound; BEAM2's # Saint-Venant kernel matches the closed form to machine precision once # the user supplies the correct K (β·b·t^3) for the cross-section. assert abs(err_pct) < 2.0, f"tip twist deviation {err_pct:.4f}% exceeds 2 %" .. rst-class:: sphx-glr-script-out .. code-block:: none φ_tip computed = 1.475177e-04 rad φ_tip published = 1.479374e-04 rad relative error = -0.2837 % .. GENERATED FROM PYTHON SOURCE LINES 105-107 Render the deformed beam ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 107-145 .. code-block:: Python dof = m.dof_map() disp = np.asarray(res.displacement, dtype=np.float64) pts = np.asarray(m.grid.points) rotx = np.zeros(pts.shape[0]) for i in range(pts.shape[0]): rows = np.where(dof[:, 0] == i + 1)[0] for r in rows: if dof[r, 1] == 3: # ROTX rotx[i] = disp[r] # Visualise the twist by superimposing a "fibre line" at one corner of # the cross-section that rotates by rotx about the beam axis as we # walk down the span. The fibre starts at (x, b/2, 0) and rotates # by phi(x) = rotx[i] about the x-axis. half_b = problem.width / 2.0 fibre_pts = np.zeros((pts.shape[0], 3)) for i, x in enumerate(pts[:, 0]): a = rotx[i] fibre_pts[i] = [x, half_b * np.cos(a), half_b * np.sin(a)] fibre = pv.PolyData(fibre_pts) # Connect consecutive points into a polyline. n = pts.shape[0] fibre.lines = np.concatenate([[n], np.arange(n)]) plotter = pv.Plotter(off_screen=True, window_size=(720, 360)) plotter.add_mesh(m.grid, color="grey", line_width=2, opacity=0.5, label="beam axis") plotter.add_mesh(fibre, color="#d62728", line_width=4, label="corner fibre (deformed)") straight = np.column_stack([pts[:, 0], np.full_like(pts[:, 0], half_b), np.zeros_like(pts[:, 0])]) straight_pd = pv.PolyData(straight) straight_pd.lines = np.concatenate([[n], np.arange(n)]) plotter.add_mesh( straight_pd, color="grey", line_width=2, opacity=0.5, label="corner fibre (undeformed)" ) plotter.add_legend() plotter.camera_position = [(2.0, -1.0, 1.0), (0.5, 0.0, 0.0), (0.0, 0.0, 1.0)] plotter.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_cantilever_torsion_001.png :alt: example verify cantilever torsion :srcset: /gallery/verification/images/sphx_glr_example_verify_cantilever_torsion_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.166 seconds) .. _sphx_glr_download_gallery_verification_example_verify_cantilever_torsion.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_torsion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_cantilever_torsion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_cantilever_torsion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_