{ "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# Pinched ring \u2014 Castigliano diametrical deflection\n\nClassical thin-circular-ring elasticity benchmark: a ring of\nmean radius $R$ and rectangular cross-section\n$b \\times t$ (with $b \\ll R$ and $t \\ll R$)\nis loaded by two equal opposing point loads $P$ along a\ndiameter. Castigliano's theorem applied to the curved-beam\nstrain energy gives the **diametrical deflection** along the\nloading axis (Timoshenko & Young 1968 \u00a779; Roark Table 9 case\n1):\n\n\\begin{align}:label: pinched-ring-castigliano\n\n \\delta_{\\mathrm{diam}}\n \\;=\\;\n \\frac{P\\, R^{3}}{E\\, I}\\,\n \\left(\\frac{\\pi}{4} - \\frac{2}{\\pi}\\right),\\end{align}\n\nwith $I = b\\, t^{3} / 12$ the curved-beam bending\ninertia about the cross-section's neutral axis perpendicular\nto the ring plane.\n\nA quarter-symmetry FE model exploits the two orthogonal\nmirror planes of the loaded ring; the loaded point's\nin-plane displacement equals $\\delta_{\\mathrm{diam}} /\n2$ since each half of the ring contributes one half of the\ndiametrical motion.\n\n## Implementation\n\nDrives the existing\n:class:`~femorph_solver.validation.problems.PinchedRing`\nproblem class. Quarter-arc of a thin ring discretised with\n40 BEAM2 elements; the loaded end has a single\n$P$-direction force; the symmetry-plane end has the\nmirror-image rotation / out-of-plane DOFs pinned.\n\nThe classical Castigliano derivation uses inextensible-ring\napproximations (no axial stretching contribution) \u2014 at the\ndefault geometry $R = 0.1$, $b \\times t = 0.01\n\\times 0.005$, the FE solution recovers the closed form to\nwithin 0.1 %, well below the 5 % tolerance the validation\nsuite specifies.\n\n## References\n\n* Timoshenko, S. P. and Young, D. H. (1968) *Elements of\n Strength of Materials*, 5th ed., Van Nostrand, \u00a779\n (Castigliano on a circular ring).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n Stress and Strain*, 6th ed., McGraw-Hill, Table 9 case 1\n (closed circular ring under two opposing point loads).\n* Bickford, W. B. (1968) *Advanced Mechanics of Materials*,\n Addison-Wesley, \u00a76.4 (curved-beam strain energy).\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.\n (2002) *Concepts and Applications of Finite Element\n Analysis*, 4th ed., Wiley, \u00a716.5 (curved-beam FE).\n\n## Vendor cross-references\n\n====================================== ========================= ============================================\nSource Reported \u03b4_diam/2 [m] Problem ID / location\n====================================== ========================= ============================================\nClosed form (Castigliano) 3.57 \u00d7 10\u207b\u2075 Timoshenko & Young 1968 \u00a779\nRoark Formulas for Stress and Strain 3.57 \u00d7 10\u207b\u2075 8th ed. Table 9.2 Case 13\nMAPDL Verification Manual \u2248 3.57 \u00d7 10\u207b\u2075 VM-38 pinched-ring family\nAbaqus Verification Manual \u2248 3.57 \u00d7 10\u207b\u2075 AVM 1.5.x pinched-ring family\nNAFEMS R0011 \u2248 3.57 \u00d7 10\u207b\u2075 pinched-ring benchmark (quarter-symmetry)\n====================================== ========================= ============================================\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport math\n\nimport numpy as np\nimport pyvista as pv\n\nfrom femorph_solver.validation.problems import PinchedRing" ] }, { "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 = PinchedRing()\nm = problem.build_model()\n\n# Quarter-arc geometry: 40 BEAM2 elements along \u03c0/2 of the ring.\nprint(\n f\"pinched ring quarter-arc: {m.grid.n_points} nodes, {m.grid.n_cells} BEAM2 cells; \"\n f\"R = {problem.R} m, b \u00d7 t = {problem.width} \u00d7 {problem.height} m, P = {problem.P} N\"\n)\nprint(f\"E = {problem.E / 1e9:.0f} GPa, \u03bd = {problem.nu}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Closed-form Castigliano reference\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "I_section = problem.width * problem.height**3 / 12.0\ndelta_full = problem.P * problem.R**3 / (problem.E * I_section) * (math.pi / 4.0 - 2.0 / math.pi)\nux_quarter = delta_full / 2.0 # quarter-symmetry split\n\nprint(f\"I (b\u00b7t\u00b3/12) = {I_section:.3e} m^4\")\nprint(f\"\u03b4_diam (full ring, Timoshenko \u00a779) = {delta_full * 1e6:.3f} \u00b5m\")\nprint(f\"UX at loaded point (quarter model) = {ux_quarter * 1e6:.3f} \u00b5m\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Static solve + extract\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res = m.solve_static()\nux_computed = problem.extract(m, res, \"loaded_point_ux\")\nerr_pct = (ux_computed - ux_quarter) / ux_quarter * 100.0\n\nprint()\nprint(f\"UX computed = {ux_computed * 1e6:8.4f} \u00b5m\")\nprint(f\"UX published = {ux_quarter * 1e6:8.4f} \u00b5m (Castigliano)\")\nprint(f\"relative error = {err_pct:+.4f} %\")\n\n# 0.5% tolerance \u2014 BEAM2 with 40 segments captures the curved-beam\n# strain energy to within ~0.1 % on this geometry; the validation\n# class permits 5 % to soak up coarser-mesh runs.\nassert abs(err_pct) < 0.5, f\"deflection deviation {err_pct:.4f}% exceeds 0.5 %\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the deformed quarter-arc + the analytical diametral motion\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)\npts = np.asarray(m.grid.points)\ndisp_xyz = np.zeros_like(pts)\nfor i in range(pts.shape[0]):\n rows = np.where(dof[:, 0] == i + 1)[0]\n for r in rows:\n d_idx = int(dof[r, 1])\n if d_idx < 3:\n disp_xyz[i, d_idx] = float(disp[r])\n\nwarped = m.grid.copy()\nwarped.points = pts + 1.0e3 * disp_xyz # \u00d71000 warp factor\n\nplotter = pv.Plotter(off_screen=True, window_size=(640, 640))\nplotter.add_mesh(m.grid, color=\"grey\", line_width=2, opacity=0.4, label=\"undeformed quarter\")\nplotter.add_mesh(warped, color=\"#1f77b4\", line_width=4, label=\"deformed (\u00d71000 warp)\")\n\n# Mark the two characteristic loci: loaded point (P applied) at\n# (R, 0) and the symmetry-plane apex (0, R).\nplotter.add_points(\n np.array([[problem.R, 0.0, 0.0]]),\n render_points_as_spheres=True,\n point_size=18,\n color=\"#d62728\",\n label=f\"loaded point \u2014 UX = {ux_computed * 1e6:.3f} \u00b5m\",\n)\nplotter.add_points(\n np.array([[0.0, problem.R, 0.0]]),\n render_points_as_spheres=True,\n point_size=14,\n color=\"black\",\n label=\"symmetry plane (UY-mirror)\",\n)\nplotter.add_legend()\nplotter.view_xy()\nplotter.camera.zoom(1.05)\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 }