{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Cantilever Saint-Venant torsion \u2014 rectangular cross-section\n\nClassical Saint-Venant elasticity verification: a clamped-free\nstraight beam carrying an end torque $T$ at its free tip\ndevelops a uniform twist rate. The total tip rotation is\n(Roark & Young 1989 Table 20 case 4; Timoshenko & Goodier 1970\n\u00a790):\n\n\\begin{align}\\varphi_\\mathrm{tip} = \\frac{T\\, L}{G\\, K},\n \\qquad\n G = \\frac{E}{2 (1 + \\nu)},\\end{align}\n\nwith $K$ the Saint-Venant torsion constant. For a\n**rectangular** cross-section $b \\times t$ (with\n$b \\ge t$) the converged Saint-Venant series gives\n\n\\begin{align}K = \\beta\\bigl(b / t\\bigr)\\, b\\, t^{3},\n \\qquad\n \\beta(1) = 0.1406,\\end{align}\n\n(Roark \u00a710 Table 20). Note: the polar moment $J = b t\n(b^{2} + t^{2}) / 12$ is **not** the torsion constant for non-\ncircular sections \u2014 using $J$ instead of $K$\noverstates the stiffness of a square shaft by ~12 %.\n\n## Implementation\n\nDrives the existing\n:class:`~femorph_solver.validation.problems.CantileverTorsion`\nproblem class on a 20-element BEAM2 (Cook Table 16.3-1 Hermite-\ncubic + Saint-Venant torsion) line. The validation suite uses\na 2 % engineering tolerance \u2014 the tip twist matches the closed\nform to machine precision.\n\n## References\n\n* Saint-Venant, B. (1855) \"M\u00e9moire sur la torsion des prismes\",\n *M\u00e9m. Acad. Sci.* 14, 233\u2013560.\n* Timoshenko, S. P. and Goodier, J. N. (1970) *Theory of\n Elasticity*, 3rd ed., McGraw-Hill, \u00a790 (rectangular\n Saint-Venant torsion).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n Stress and Strain*, 6th ed., McGraw-Hill, Table 20 case 4\n (rectangular Saint-Venant torsion constant).\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n *Concepts and Applications of Finite Element Analysis*, 4th\n ed., Wiley, Table 16.3-1 (BEAM2 Saint-Venant kernel).\n\n## Vendor cross-references\n\n================================ ===================== ===================================\nSource Reported \u03c6_tip [rad] Problem ID / location\n================================ ===================== ===================================\nClosed form (Saint-Venant) 1.475 \u00d7 10\u207b\u2074 Timoshenko & Goodier 1970 \u00a7109\nSaint-Venant (1853) 1.475 \u00d7 10\u207b\u2074 original derivation, \u03b2(b/t) series\nMAPDL Verification Manual 1.475 \u00d7 10\u207b\u2074 VM-23 family (torsion kernel)\nAbaqus Verification Manual 1.475 \u00d7 10\u207b\u2074 AVM 1.5.x cantilever-torsion family\n================================ ===================== ===================================\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport numpy as np\nimport pyvista as pv\n\nfrom femorph_solver.validation.problems import CantileverTorsion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build the model from the validation problem class\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "problem = CantileverTorsion()\nm = problem.build_model()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Closed-form quantities\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "beta = 0.1406 # square cross-section (b/t = 1)\nG = problem.E / (2.0 * (1.0 + problem.nu))\nK = beta * problem.width * problem.height**3\nphi_published = problem.T * problem.L / (G * K)\n\nprint(f\"L = {problem.L} m, b \u00d7 t = {problem.width} \u00d7 {problem.height} m, T = {problem.T} N m\")\nprint(f\"E = {problem.E / 1e9:.0f} GPa, \u03bd = {problem.nu}, G = {G / 1e9:.2f} GPa\")\nprint(f\"\u03b2(b/t = 1) = {beta} (Roark Table 20 case 4)\")\nprint(f\"K = \u03b2\u00b7b\u00b7t^3 = {K:.4e} m^4 (Saint-Venant torsion constant)\")\nprint(f\"\u03c6_tip = T L / (G K) = {phi_published:.6e} rad ({np.degrees(phi_published):.4f} \u00b0)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Static solve + tip-twist extraction\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res = m.solve_static()\nphi_computed = problem.extract(m, res, \"tip_twist\")\nerr_pct = (phi_computed - phi_published) / phi_published * 100.0\nprint()\nprint(f\"\u03c6_tip computed = {phi_computed:.6e} rad\")\nprint(f\"\u03c6_tip published = {phi_published:.6e} rad\")\nprint(f\"relative error = {err_pct:+.4f} %\")\n\n# 2% tolerance \u2014 the validation class's tight engineering bound; BEAM2's\n# Saint-Venant kernel matches the closed form to machine precision once\n# the user supplies the correct K (\u03b2\u00b7b\u00b7t^3) for the cross-section.\nassert abs(err_pct) < 2.0, f\"tip twist deviation {err_pct:.4f}% exceeds 2 %\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the deformed beam\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dof = m.dof_map()\ndisp = np.asarray(res.displacement, dtype=np.float64)\n\npts = np.asarray(m.grid.points)\nrotx = np.zeros(pts.shape[0])\nfor i in range(pts.shape[0]):\n rows = np.where(dof[:, 0] == i + 1)[0]\n for r in rows:\n if dof[r, 1] == 3: # ROTX\n rotx[i] = disp[r]\n\n# Visualise the twist by superimposing a \"fibre line\" at one corner of\n# the cross-section that rotates by rotx about the beam axis as we\n# walk down the span. The fibre starts at (x, b/2, 0) and rotates\n# by phi(x) = rotx[i] about the x-axis.\nhalf_b = problem.width / 2.0\nfibre_pts = np.zeros((pts.shape[0], 3))\nfor i, x in enumerate(pts[:, 0]):\n a = rotx[i]\n fibre_pts[i] = [x, half_b * np.cos(a), half_b * np.sin(a)]\nfibre = pv.PolyData(fibre_pts)\n# Connect consecutive points into a polyline.\nn = pts.shape[0]\nfibre.lines = np.concatenate([[n], np.arange(n)])\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 360))\nplotter.add_mesh(m.grid, color=\"grey\", line_width=2, opacity=0.5, label=\"beam axis\")\nplotter.add_mesh(fibre, color=\"#d62728\", line_width=4, label=\"corner fibre (deformed)\")\nstraight = np.column_stack([pts[:, 0], np.full_like(pts[:, 0], half_b), np.zeros_like(pts[:, 0])])\nstraight_pd = pv.PolyData(straight)\nstraight_pd.lines = np.concatenate([[n], np.arange(n)])\nplotter.add_mesh(\n straight_pd, color=\"grey\", line_width=2, opacity=0.5, label=\"corner fibre (undeformed)\"\n)\nplotter.add_legend()\nplotter.camera_position = [(2.0, -1.0, 1.0), (0.5, 0.0, 0.0), (0.0, 0.0, 1.0)]\nplotter.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 0 }