{ "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# NAFEMS LE1 \u2014 elliptic membrane (plane stress)\n\nCanonical plane-stress benchmark from the NAFEMS *Standard\nBenchmark Tests for Linear Elastic Analysis*, R0015 (1990),\nproblem **LE1**. A quarter elliptic-membrane plate bounded by\ntwo confocal ellipses\n\n\\begin{align}\\frac{x^{2}}{2^{2}} + \\frac{y^{2}}{1^{2}} = 1\n \\qquad \\text{(inner)},\\end{align}\n\n\\begin{align}\\frac{x^{2}}{3.25^{2}} + \\frac{y^{2}}{2.75^{2}} = 1\n \\qquad \\text{(outer)},\\end{align}\n\nis loaded by a uniform 1 MPa outward-normal pressure on the\nouter elliptic edge. Symmetry pins $u_y = 0$ on the\n$x$-axis and $u_x = 0$ on the $y$-axis.\n\nThe published reference is the **hoop stress at point D**\n$(x = 2,\\ y = 0)$ \u2014 the stub of the inner ellipse on the\n$x$-axis:\n\n\\begin{align}\\sigma_{yy}^{(D)} = 92.7 \\;\\text{MPa}.\\end{align}\n\nNo analytical closed form exists; the value is benchmark-defined\nfrom converged solutions reported by multiple independent codes\nin the original 1988 NAFEMS working-group validation.\n\n## Implementation\n\nThis gallery example drives the same\n:class:`~femorph_solver.validation.problems.NafemsLE1` problem\nclass that the validation suite exercises, on a moderate\n$16 \\times 4 = 64$-cell QUAD4 plane-stress mesh. The\n\u03c3_yy(D) recovery uses a finite-difference of the inner-ellipse\nnodal displacement field, in lock-step with the workaround the\nvalidation problem class itself ships.\n\n## References\n\n* NAFEMS, *Standard Benchmark Tests for Linear Elastic\n Analysis*, R0015 (1990), \u00a72.1 LE1.\n* NAFEMS, *Background to Benchmarks*, R0013 (1988) \u2014\n methodology companion describing the converged reference.\n* NAFEMS LE1 plane-stress membrane (canonical reference); Abaqus AVM\n 1.3.x elliptic-membrane plane-stress family.\n\n## Vendor cross-references\n\n====================================== ========================= ============================================\nSource Reported \u03c3_yy at D [MPa] Problem ID / location\n====================================== ========================= ============================================\nNAFEMS R0015 \u00a72.1 92.7 LE1 Elliptic membrane under outward pressure\nNAFEMS R0013 (methodology) 92.7 Background to Benchmarks companion\nMAPDL Verification Manual \u2248 92.7 NAFEMS LE1 plane-stress membrane entry\nAbaqus Verification Manual \u2248 92.7 AVM 1.3.x elliptic-membrane 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 NafemsLE1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Build the model from the validation problem class\n\nParameters (semi-axes / pressure / E / \u03bd) come from\n:class:`NafemsLE1`'s defaults; the problem class assembles a\nQUAD4_PLANE mapped mesh and applies the consistent outer-edge\npressure ring.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "problem = NafemsLE1()\nm = problem.build_model(n_theta=16, n_rad=4)\nprint(\n f\"NAFEMS LE1 mesh: {m.grid.n_points} nodes, {m.grid.n_cells} QUAD4 cells; \"\n f\"thickness = {problem.thickness} m, p_outer = {problem.pressure / 1e6:.1f} MPa\"\n)\nprint(f\"E = {problem.E / 1e9:.0f} GPa, \u03bd = {problem.nu}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Static solve\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res = m.solve_static()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## \u03c3_yy at point D \u2014 finite-difference recovery\n\nPoint D is the inner-ellipse node at $(2, 0)$. Plane-\nstress constitutive law gives\n\n\\begin{align}\\sigma_{yy} = \\frac{E}{1 - \\nu^{2}}\n (\\varepsilon_{yy} + \\nu\\, \\varepsilon_{xx}),\\end{align}\n\nwith strains estimated by forward finite-differences against\nthe next radial / \u03b8 nodes on the parametric grid.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sigma_yy_d = problem.extract(m, res, \"sigma_yy_at_D\")\nsigma_yy_d_published = 92.7e6\nerr_pct = (sigma_yy_d - sigma_yy_d_published) / sigma_yy_d_published * 100.0\n\nprint()\nprint(f\"\u03c3_yy at D computed = {sigma_yy_d / 1e6:.3f} MPa\")\nprint(f\"\u03c3_yy at D published = {sigma_yy_d_published / 1e6:.3f} MPa (NAFEMS R0015 LE1)\")\nprint(f\"relative error = {err_pct:+.2f} %\")\n\n# Allow the same 8 % tolerance the validation problem class uses.\nassert abs(err_pct) < 8.0, f\"\u03c3_yy(D) deviation {err_pct:.2f}% exceeds NAFEMS engineering tolerance\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the deformed mesh, coloured by displacement magnitude\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid = m.grid.copy()\nu = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 2)\ndisp3 = np.zeros((grid.n_points, 3), dtype=np.float64)\ndisp3[:, :2] = u\ngrid.point_data[\"displacement\"] = disp3\ngrid.point_data[\"|u|\"] = np.linalg.norm(u, axis=1)\n\nwarped = grid.warp_by_vector(\"displacement\", factor=2.0e3)\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 480))\nplotter.add_mesh(\n grid,\n style=\"wireframe\",\n color=\"grey\",\n opacity=0.4,\n line_width=1,\n label=\"undeformed\",\n)\nplotter.add_mesh(\n warped,\n scalars=\"|u|\",\n cmap=\"plasma\",\n show_edges=True,\n scalar_bar_args={\"title\": \"|u| [m] (deformed \u00d72000)\"},\n)\n\n# Mark point D and the y-axis (symmetry plane) for orientation.\nplotter.add_points(\n np.array([[problem.a_in, 0.0, 0.0]]),\n render_points_as_spheres=True,\n point_size=18,\n color=\"#d62728\",\n label=\"point D = (2, 0)\",\n)\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 }