{
  "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# Tutorial 3 \u2014 Pressure-vessel design-by-analysis (Lam\u00e9 thick cylinder)\n\nYou're sizing the wall of a small pressure vessel that operates\nat $p = 10\\,\\text{MPa}$ internal pressure.  The\npreliminary design has bore radius $a = 100\\,\\text{mm}$\nand outer radius $b = 200\\,\\text{mm}$ \u2014 a thick-walled\ngeometry where the simpler thin-wall hoop-stress approximation\n$\\sigma_\\theta = p\\, r / t$ doesn't apply.  ASME B31.3\n\u00a7304.1 mandates that for $r/t < 6$ you have to use either\nthe full Lam\u00e9 closed form or a numerical analysis.\n\nThis tutorial walks the design-by-analysis loop end-to-end on\nthe Lam\u00e9 thick cylinder.  Steps:\n\n* **Step 1** \u2014 build a quarter-symmetry HEX8 model of the wall.\n* **Step 2** \u2014 apply the plane-strain symmetry BCs and the\n  internal-pressure load distribution.\n* **Step 3** \u2014 static solve.\n* **Step 4** \u2014 recover the stress field; tabulate\n  $\\sigma_\\theta$ and $\\sigma_r$ at the bore, mid-\n  radius, and outer surface against the Lam\u00e9 closed form.\n* **Step 5** \u2014 compute $\\sigma_\\text{VM}$ at the bore and\n  compare against the ASME design allowable.\n* **Step 6** \u2014 refine the mesh and confirm the bore-edge stress\n  has converged.\n* **Step 7** \u2014 render the deformed pressure vessel coloured by\n  $\\sigma_\\text{VM}$.\n\nThe companion focused recipes are\n`sphx_glr_gallery_post-processing_example_nodal_stress_recovery.py`\n(stress recovery + invariants on the same geometry, narrower\nscope) and\n`sphx_glr_gallery_verification_example_verify_lame_cylinder.py`\n(pure verification \u2014 same closed form, no design check).\n\n## Theoretical reference\n\nFor a thick cylinder under internal pressure $p_i$ and\nzero external pressure, Lam\u00e9's 1852 closed form\n(Timoshenko & Goodier 1970 \u00a733) gives the radial and hoop\nstresses\n\n\\begin{align}\\sigma_r(r) = -\\frac{p_i\\, a^2\\, (b^2 - r^2)}{r^2\\, (b^2 - a^2)},\n    \\qquad\n    \\sigma_\\theta(r) = +\\frac{p_i\\, a^2\\, (b^2 + r^2)}{r^2\\, (b^2 - a^2)}.\\end{align}\n\nAt the bore ($r = a$), $\\sigma_r(a) = -p_i$ (free-\nsurface boundary condition) and the hoop stress peaks at\n\n\\begin{align}\\sigma_\\theta(a)\n        = p_i\\, \\frac{a^2 + b^2}{b^2 - a^2}.\\end{align}\n\nFor our $a = 0.1\\,\\text{m}, b = 0.2\\,\\text{m},\np_i = 10\\,\\text{MPa}$ design, the closed form gives\n$\\sigma_\\theta(a) \\approx 16.67\\,\\text{MPa}$ \u2014\nsubstantially below the thin-wall guess $pa/t = 10\\,\n\\text{MPa}$.  Lam\u00e9 is *less* conservative than thin-wall here.\n\nDesign-allowable convention: ASME B31.3 \u00a7302.3 specifies the\nbasic allowable stress $S$ at design temperature.  For\nSA-516 Gr. 70 carbon steel at room temperature,\n$S \\approx 138\\,\\text{MPa}$; we'll use $S$ as the\nallowable $\\sigma_\\text{VM}$.\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\nimport femorph_solver\nfrom femorph_solver import ELEMENTS\nfrom femorph_solver.recover import compute_nodal_stress, stress_invariants"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 1 \u2014 build a quarter-symmetry HEX8 model of the wall\n\nA pressure vessel with no caps and no body-force gradient is\naxisymmetric, so a quarter-symmetry model captures the full\nresponse.  We build the wall as a structured \u03b8-r-z grid then\nconvert to HEX8 cells.  The slab thickness $t_z = 20\\,\n\\text{mm}$ is held fully in the plane-strain sense by the BCs\nin Step 2.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = 0.10  # bore radius [m]\nb = 0.20  # outer radius [m]\nt_axial = 0.02  # plane-strain slab thickness [m]\np_i = 10.0e6  # internal pressure [Pa]\nE = 2.0e11  # SA-516 Gr 70 steel\nNU = 0.30\nRHO = 7850.0\nsigma_allow = 138.0e6  # ASME B31.3 \u00a7302.3 basic allowable for SA-516 Gr 70 at RT\n\n# Closed-form predictions used for comparison\nsigma_theta_a_pub = p_i * (a**2 + b**2) / (b**2 - a**2)\nsigma_r_a_pub = -p_i\n\nprint(\"Pressure-vessel design-by-analysis \u2014 Lam\u00e9 thick cylinder\")\nprint(f\"  a = {a} m, b = {b} m  \u2192  thick-wall ratio b/a = {b / a:.2f}\")\nprint(f\"  p_i = {p_i / 1e6:.1f} MPa, E = {E / 1e9:.0f} GPa, \u03bd = {NU}\")\nprint(f\"  \u03c3_allow (ASME B31.3, SA-516 Gr 70 RT) = {sigma_allow / 1e6:.0f} MPa\")\nprint()\nprint(\"  Closed form (Lam\u00e9 1852):\")\nprint(f\"    \u03c3_\u03b8(a) = p_i\u00b7(a\u00b2+b\u00b2)/(b\u00b2-a\u00b2) = {sigma_theta_a_pub / 1e6:.3f} MPa\")\nprint(f\"    \u03c3_r(a) = -p_i               = {sigma_r_a_pub / 1e6:.3f} MPa\")\n\n\ndef build_wall(n_theta: int, n_rad: int) -> tuple[femorph_solver.Model, np.ndarray]:\n    \"\"\"Build a quarter-annulus HEX8-EAS plane-strain mesh.\"\"\"\n    theta = np.linspace(0.0, 0.5 * math.pi, n_theta + 1)\n    r = np.linspace(a, b, n_rad + 1)\n\n    pts_list: list[list[float]] = []\n    for kz in (0.0, t_axial):\n        for ti in theta:\n            for rj in r:\n                pts_list.append([rj * math.cos(ti), rj * math.sin(ti), kz])\n    pts_arr = np.array(pts_list, dtype=np.float64)\n\n    nx_plane = (n_theta + 1) * (n_rad + 1)\n    n_cells = n_theta * n_rad\n    cells = np.empty((n_cells, 9), dtype=np.int64)\n    cells[:, 0] = 8\n    c = 0\n    for i in range(n_theta):\n        for j in range(n_rad):\n            n00 = i * (n_rad + 1) + j\n            n10 = i * (n_rad + 1) + (j + 1)\n            n11 = (i + 1) * (n_rad + 1) + (j + 1)\n            n01 = (i + 1) * (n_rad + 1) + j\n            cells[c, 1:] = (\n                n00,\n                n10,\n                n11,\n                n01,\n                n00 + nx_plane,\n                n10 + nx_plane,\n                n11 + nx_plane,\n                n01 + nx_plane,\n            )\n            c += 1\n    grid = pv.UnstructuredGrid(cells.ravel(), np.full(n_cells, 12, dtype=np.uint8), pts_arr)\n\n    m = femorph_solver.Model.from_grid(grid)\n    m.assign(\n        ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n        material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n    )\n    return m, pts_arr"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 2 \u2014 boundary conditions and pressure load\n\n*Plane-strain* + quarter-symmetry needs three BC patches:\n\n- On the $\\theta = 0$ face ($y = 0$), pin\n  $u_y = 0$ \u2014 the symmetry condition.\n- On the $\\theta = \\pi/2$ face ($x = 0$), pin\n  $u_x = 0$.\n- On every node, pin $u_z = 0$ to enforce plane strain\n  (no axial displacement).\n\nThe internal-pressure load is **lumped onto the bore-face\nnodes** by integrating the pressure over each \u03b8-segment via\nthe trapezoid rule.  Per-segment force = $p_i \\times ds\n\\times t_\\text{slab}$, distributed half to each segment\nendpoint along the outward radial direction.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def apply_bcs_and_load(m: femorph_solver.Model, pts_arr: np.ndarray, n_theta: int, n_rad: int):\n    nx_plane = (n_theta + 1) * (n_rad + 1)\n    for k, p in enumerate(pts_arr):\n        if p[0] < 1e-9:  # x = 0 face\n            m.fix(nodes=[int(k + 1)], dof=\"UX\")\n        if p[1] < 1e-9:  # y = 0 face\n            m.fix(nodes=[int(k + 1)], dof=\"UY\")\n        m.fix(nodes=[int(k + 1)], dof=\"UZ\")  # plane strain everywhere\n\n    # Internal-pressure load at the bore (j = 0 in the radial\n    # ladder).  Loop over \u03b8-segments and lump the per-segment\n    # pressure resultant onto the two endpoints.\n    fx_acc: dict[int, float] = {}\n    fy_acc: dict[int, float] = {}\n    for kz in (0, 1):\n        inner = [kz * nx_plane + i * (n_rad + 1) + 0 for i in range(n_theta + 1)]\n        for seg in range(n_theta):\n            ai, bi = inner[seg], inner[seg + 1]\n            ds = float(np.linalg.norm(pts_arr[bi] - pts_arr[ai]))\n            mid = 0.5 * (pts_arr[ai] + pts_arr[bi])\n            rxy = np.array([mid[0], mid[1]])\n            nrm = float(np.linalg.norm(rxy))\n            outward = rxy / nrm if nrm > 1e-12 else np.zeros(2)\n            f_seg = p_i * ds * (t_axial / 2.0)\n            for n_idx in (ai, bi):\n                fx_acc[n_idx] = fx_acc.get(n_idx, 0.0) + 0.5 * f_seg * outward[0]\n                fy_acc[n_idx] = fy_acc.get(n_idx, 0.0) + 0.5 * f_seg * outward[1]\n    for n_idx, fx in fx_acc.items():\n        fy = fy_acc[n_idx]\n        m.apply_force(int(n_idx + 1), fx=fx, fy=fy)\n\n\nN_THETA, N_RAD = 24, 16\nm, pts_arr = build_wall(N_THETA, N_RAD)\napply_bcs_and_load(m, pts_arr, N_THETA, N_RAD)\nprint(f\"\\n  Mesh: ({N_THETA} \u00d7 {N_RAD}) HEX8-EAS = {m.grid.n_cells} cells, {m.grid.n_points} nodes\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 3 \u2014 static solve\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\nprint(f\"  {res!r}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 4 \u2014 stress recovery and Lam\u00e9 tabulation\n\nExtract the per-node Voigt stress and the principal-stress\nscalars.  At points along the $y = 0$ line, the radial\ndirection coincides with $+x$ and the hoop direction\nwith $+y$, so $\\sigma_r$ reads as $\\sigma_{xx}$\nand $\\sigma_\\theta$ reads as $\\sigma_{yy}$ \u2014 handy\nfor tabulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "u = np.asarray(res.displacement, dtype=np.float64).ravel()\nsigma = compute_nodal_stress(m, u)\ninv = stress_invariants(sigma)\n\n\ndef lookup_radial(target_r: float) -> tuple[float, float, float]:\n    \"\"\"Return (\u03c3_xx, \u03c3_yy, \u03c3_VM) at the equator point closest to r = target.\"\"\"\n    diff = np.linalg.norm(pts_arr[:, :2] - np.array([target_r, 0.0]), axis=1)\n    i_star = int(np.argmin(diff))\n    return (\n        float(sigma[i_star, 0]),\n        float(sigma[i_star, 1]),\n        float(inv[\"von_mises\"][i_star]),\n    )\n\n\nprint()\nprint(\n    f\"  {'r [m]':>7}  {'\u03c3_r FE [MPa]':>14}  {'\u03c3_r pub [MPa]':>14}  \"\n    f\"{'\u03c3_\u03b8 FE [MPa]':>14}  {'\u03c3_\u03b8 pub [MPa]':>14}  {'\u03c3_VM [MPa]':>11}\"\n)\nprint(\"  \" + \"-\" * 90)\nfor r_target in (a, 0.5 * (a + b), b):\n    sxx, syy, svm = lookup_radial(r_target)\n    sigma_theta_pub = p_i * a**2 * (b**2 + r_target**2) / (r_target**2 * (b**2 - a**2))\n    sigma_r_pub = -p_i * a**2 * (b**2 - r_target**2) / (r_target**2 * (b**2 - a**2))\n    print(\n        f\"  {r_target:>7.4f}  {sxx / 1e6:>+12.3f}  {sigma_r_pub / 1e6:>+12.3f}  \"\n        f\"{syy / 1e6:>+12.3f}  {sigma_theta_pub / 1e6:>+12.3f}  {svm / 1e6:>10.3f}\"\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 5 \u2014 engineering safety check\n\nThe peak von-Mises stress lives at the bore ($r = a$) by\ninspection.  Compare against the ASME basic allowable\n$S$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "_, _, sigma_vm_at_bore = lookup_radial(a)\nsafety_factor = sigma_allow / sigma_vm_at_bore\nprint()\nprint(f\"  Peak \u03c3_VM at bore = {sigma_vm_at_bore / 1e6:.3f} MPa\")\nprint(f\"  \u03c3_allow (ASME B31.3) = {sigma_allow / 1e6:.0f} MPa\")\nprint(f\"  Safety factor = \u03c3_allow / \u03c3_VM_max = {safety_factor:.2f}\")\nif safety_factor >= 1.5:\n    print(f\"  PASS: SF = {safety_factor:.2f} \u2265 1.5 design margin.\")\nelse:\n    print(f\"  FAIL: SF = {safety_factor:.2f} < 1.5 \u2014 increase wall thickness.\")\nassert safety_factor >= 1.5, \"design fails the 1.5 SF check\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 6 \u2014 mesh-refinement convergence\n\nThe bore-edge stress is the most refinement-sensitive value\nwe recover \u2014 the \u03c3_\u03b8 gradient steepens as $r \\to a$,\nand a coarse mesh can over- or under-predict the peak.  Sweep\nthree meshes and read \u03c3_\u03b8(a) off each.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(\n    f\"  {'mesh (N_\u03b8, N_r)':>16}  {'\u03c3_\u03b8(a) FE [MPa]':>16}  {'err vs Lam\u00e9 [%]':>17}  \"\n    f\"{'\u03c3_VM peak [MPa]':>16}\"\n)\nprint(\"  \" + \"-\" * 75)\nfor n_th, n_r in ((12, 8), (24, 16), (36, 24)):\n    m_r, pts_r = build_wall(n_th, n_r)\n    apply_bcs_and_load(m_r, pts_r, n_th, n_r)\n    res_r = m_r.solve_static()\n    u_r = np.asarray(res_r.displacement, dtype=np.float64).ravel()\n    sigma_r_field = compute_nodal_stress(m_r, u_r)\n    inv_r = stress_invariants(sigma_r_field)\n    diff = np.linalg.norm(pts_r[:, :2] - np.array([a, 0.0]), axis=1)\n    i_bore = int(np.argmin(diff))\n    syy_bore = float(sigma_r_field[i_bore, 1])\n    svm_peak = float(inv_r[\"von_mises\"].max())\n    err_pct = (syy_bore - sigma_theta_a_pub) / sigma_theta_a_pub * 100.0\n    print(\n        f\"  ({n_th:>3}, {n_r:>3})    {syy_bore / 1e6:>+14.3f}     {err_pct:>+15.3f}    \"\n        f\"{svm_peak / 1e6:>15.3f}\"\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Step 7 \u2014 render \u03c3_VM on the deformed wall\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid_render = m.grid.copy()\ngrid_render.point_data[\"\u03c3_VM [Pa]\"] = inv[\"von_mises\"]\ngrid_render.point_data[\"displacement\"] = u.reshape(-1, 3)\nwarped = grid_render.warp_by_vector(\"displacement\", factor=2.0e3)\n\nplotter = pv.Plotter(off_screen=True, window_size=(820, 540))\nplotter.add_mesh(grid_render, style=\"wireframe\", color=\"grey\", opacity=0.35, line_width=1)\nplotter.add_mesh(\n    warped,\n    scalars=\"\u03c3_VM [Pa]\",\n    cmap=\"plasma\",\n    show_edges=False,\n    scalar_bar_args={\"title\": \"\u03c3_VM [Pa]  (deformed \u00d72 000)\"},\n)\nplotter.view_xy()\nplotter.camera.zoom(1.05)\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Take-aways\n\n- **Thick-wall geometry** (b/a < 6 in our case 2.0) needs the\n  full Lam\u00e9 closed form, not the thin-wall hoop-stress\n  approximation \u2014 the FE result confirms this.\n- **Plane-strain 3-D solid** is the right discretisation when\n  the vessel is much longer than the wall thickness; the BCs\n  are $u_z = 0$ at every node plus quarter-symmetry pins\n  on the two cut faces.\n- **Stress recovery** uses the same ``compute_nodal_stress`` +\n  ``stress_invariants`` pipeline that the post-processing\n  recipes demonstrate; here we apply it to a closed-form\n  benchmark so every recovered value can be cross-checked.\n- **The ASME safety factor** is a single ratio against the\n  basic allowable $S$ from the code's stress table.\n- **Mesh-refinement sensitivity** is highest at the bore where\n  $\\sigma_\\theta$ peaks; two-fold refinement drops the\n  bore error from ~3.5 % to ~1.5 %.  A converged answer is\n  the pre-condition for any defensible engineering report.\n\nNext steps:\n\n- Tutorial 4 (planned) \u2014 random-vibration / DDAM analysis on a\n  marine-shock-loaded equipment skid.\n- Tutorial 5 (planned) \u2014 cyclic-symmetry rotor modal survey.\n\n"
      ]
    }
  ],
  "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
}