{ "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# Mesh-refinement convergence \u2014 Lam\u00e9 thick cylinder\n\nCompanion to the\n`sphx_glr_gallery_verification_example_verify_convergence_cantilever_eb.py`\nbeam study, this example shows the same $O(h^{2})$\ndisplacement-norm convergence rate on a **3D axisymmetric\nelasticity** problem: the Lam\u00e9 thick cylinder under internal\npressure. Different physics (purely radial deformation, no\nbending), different element geometry (annular sector instead\nof beam-aspect prism), but the same Galerkin error theorem\nholds \u2014 Strang & Fix (2008) \u00a73.7 promises\n$\\|u^{h} - u\\|_{L^{2}} = O(h^{p+1})$ with $p = 1$\nfor trilinear hex.\n\nThe bore radial displacement\n$u_r(a) = (p_i\\, a^{2}\\, a / [E (b^{2} - a^{2})]) \\cdot\n[(1 - \\nu - 2\\nu^{2}) + b^{2} (1 + \\nu) / a^{2}]$\n(Lam\u00e9 1852; Timoshenko & Goodier 1970 \u00a733) gives the\nunambiguous reference value across every refinement.\n\n## Implementation\n\nQuarter-annulus HEX8 EAS plane-strain model (one axial layer)\nwith a refinement ladder\n$(N_\\theta, N_r) \\in \\{(4, 2), (8, 4), (12, 8), (16, 12),\n(24, 16)\\}$. At each refinement we run the same internal-\npressure-applied static solve as the\n`sphx_glr_gallery_verification_example_verify_lame_cylinder.py`\nsingle-mesh example, then compare\n$u_r(a)$ against the closed form. Mesh size $h$\nis taken as $\\max(\\Delta\\theta\\, a, \\Delta r)$ to\ncharacterise the *largest* element edge length on the bore-\nadjacent ring (where the gradient is steepest).\n\n## References\n\n* Lam\u00e9, G. (1852) *Le\u00e7ons sur la Th\u00e9orie Math\u00e9matique de\n l'\u00c9lasticit\u00e9 des Corps Solides*, Bachelier \u2014 original\n derivation of the radial displacement.\n* Timoshenko, S. P. and Goodier, J. N. (1970) *Theory of\n Elasticity*, 3rd ed., McGraw-Hill, \u00a733 (thick-walled\n cylinder).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n Stress and Strain*, 6th ed., McGraw-Hill, Table 32 case 1b.\n* Strang, G. and Fix, G. (2008) *An Analysis of the Finite\n Element Method*, 2nd ed., Wellesley-Cambridge, \u00a73.7\n (a-priori error estimates).\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, \u00a76 (convergence).\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed\n assumed strain methods \u2026\" *International Journal for\n Numerical Methods in Engineering* 29 (8), 1595\u20131638.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport math\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pyvista as pv\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem data\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\np_i = 1.0e7\nE = 2.0e11\nNU = 0.30\nRHO = 7850.0\n\nur_a_published = (\n p_i * a**2 * a / (E * (b**2 - a**2)) * ((1 - NU - 2 * NU**2) + b**2 * (1 + NU) / a**2)\n)\nprint(f\"Lam\u00e9 bore u_r(a) (Timoshenko & Goodier \u00a733) = {ur_a_published * 1e6:.4f} \u00b5m\")\n\n\ndef run_one(n_theta: int, n_rad: int) -> tuple[float, float]:\n \"\"\"One quarter-annulus solve. Returns (h, |relative error|).\"\"\"\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[float]] = []\n for kz in (0.0, t_axial):\n for ti in theta:\n for rj in r:\n pts.append([rj * math.cos(ti), rj * math.sin(ti), kz])\n pts_arr = np.array(pts, 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 n00b = i * (n_rad + 1) + j\n n10b = i * (n_rad + 1) + (j + 1)\n n11b = (i + 1) * (n_rad + 1) + (j + 1)\n n01b = (i + 1) * (n_rad + 1) + j\n cells[c, 1:] = (\n n00b,\n n10b,\n n11b,\n n01b,\n n00b + nx_plane,\n n10b + nx_plane,\n n11b + nx_plane,\n n01b + nx_plane,\n )\n c += 1\n grid = pv.UnstructuredGrid(\n cells.ravel(),\n np.full(n_cells, 12, dtype=np.uint8),\n pts_arr,\n )\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\n for k, p in enumerate(pts_arr):\n if p[0] < 1e-9:\n m.fix(nodes=[int(k + 1)], dof=\"UX\")\n if p[1] < 1e-9:\n m.fix(nodes=[int(k + 1)], dof=\"UY\")\n m.fix(nodes=[int(k + 1)], dof=\"UZ\")\n\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 res = m.solve_static()\n u_flat = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 3)\n inner_x = 0 # i=0, j=0, kz=0 sits at (a, 0, 0)\n ur_a_computed = float(u_flat[inner_x, 0])\n\n # Largest bore-adjacent edge. \u0394\u03b8 \u00b7 a (arc) vs \u0394r.\n h = max((0.5 * math.pi / n_theta) * a, (b - a) / n_rad)\n rel_err = abs(ur_a_computed - ur_a_published) / abs(ur_a_published)\n return h, rel_err" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Refinement ladder\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ladder = [(4, 2), (8, 4), (12, 8), (16, 12), (24, 16)]\nhs: list[float] = []\nerrs: list[float] = []\n\nprint()\nprint(f\"{'(N_\u03b8, N_r)':>11} {'h [m]':>10} {'u_r(a) FE / pub':>16} {'rel err':>10}\")\nprint(f\"{'-' * 11:>11} {'-' * 10:>10} {'-' * 16:>16} {'-' * 10:>10}\")\nfor n_th, n_r in ladder:\n h, err = run_one(n_th, n_r)\n hs.append(h)\n errs.append(err)\n print(f\"({n_th:>3}, {n_r:>3}) {h:10.5f} {(1 - err):>16.6f} {err:10.4e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Asymptotic convergence rate\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "log_h = np.log(np.asarray(hs[1:]))\nlog_err = np.log(np.asarray(errs[1:]))\nslope, intercept = np.polyfit(log_h, log_err, 1)\np_estimated = float(slope)\n\nprint()\nprint(f\"least-squares fit (skip first): rate p \u2248 {p_estimated:.3f}\")\nprint(\"expected for HEX8 EAS displacement norm (Strang & Fix \u00a73.7): p = 2\")\n\nassert 1.4 < p_estimated, f\"convergence rate {p_estimated:.3f} unexpectedly poor (< 1.4)\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the convergence plot\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(6.4, 4.0))\nax.loglog(hs, errs, \"o-\", color=\"#2ca02c\", lw=2, label=\"HEX8 EAS \u2014 Lam\u00e9 bore u_r(a)\")\n\nh_ref = np.array([hs[0], hs[-1]])\nerr_ref_p2 = errs[-1] * (h_ref / hs[-1]) ** 2\nax.loglog(h_ref, err_ref_p2, \"--\", color=\"#d62728\", lw=1.5, label=r\"reference $\\propto h^{2}$\")\n\nax.set_xlabel(r\"max bore-adjacent edge length $h$\")\nax.set_ylabel(r\"$|u_r^{h}(a) - u_r(a)| / |u_r(a)|$\")\nax.set_title(\n f\"Lam\u00e9 thick cylinder \u2014 fitted rate p \u2248 {p_estimated:.2f}\",\n fontsize=11,\n)\nax.legend(loc=\"lower right\", fontsize=9)\nax.grid(True, which=\"both\", ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verify monotone decrease\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "assert all(e > 0 for e in errs), \"errors must be positive\"\nfor prev, nxt in zip(errs[:-1], errs[1:]):\n assert nxt < prev * 1.1, \"error grew between successive refinements\"\nprint()\nprint(\"OK \u2014 error decreases monotonically with refinement.\")" ] } ], "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 }