{ "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# Principal stresses + principal directions on a thick cylinder\n\nWhere :doc:`example_nodal_stress_recovery` answered *how big* the\nstress is, this example answers *which way* it points. The Cauchy\nstress tensor $\\sigma_{ij}$ at every node has three orthogonal\n**principal directions** along which there is no shear, with the\ncorresponding **principal stresses** $\\sigma_{1} \\ge \\sigma_{2}\n\\ge \\sigma_{3}$ along those axes. Tracing those eigenvectors over a\ndomain gives the *stress-trajectory* picture beloved by mechanical\ndesigners (Timoshenko & Goodier 1970 \u00a71.7) \u2014 the visualisation that\nshows where the load actually flows.\n\nFor a thick cylinder under internal pressure, the answer is known\nexactly: the principal directions are **radial** and **hoop**, and\n$\\sigma_{1}$ is the hoop component everywhere. So the\nrecovered eigenvectors give a clean accuracy check \u2014 they should\npoint along the local $\\hat{\\theta}$ to within mesh-\ndiscretisation error.\n\nThis example uses the public stress-recovery pipeline plus\n:func:`numpy.linalg.eigh` per node to extract the principal directions\n(``stress_invariants`` returns the principal *stresses* only):\n\n1. ``compute_nodal_stress(model, u)`` \u2014 node-averaged Voigt array.\n2. ``stress_invariants(sigma)`` \u2014 $\\sigma_{1,2,3}$ and\n companion scalars per node.\n3. ``np.linalg.eigh`` on the per-node 3x3 Cauchy tensor \u2014 yields\n sorted eigenvalues *and* eigenvectors so we can identify\n *which* principal axis aligns with which physical direction.\n\n## Verification\n\nAt every node of the Lam\u00e9 annulus, the recovered $\\sigma_{1}$\ndirection must be tangent to the local hoop circle (perpendicular to\nthe radial unit vector). We compute the angle between the recovered\n$\\sigma_{1}$ axis and the analytical hoop direction\n$\\hat{\\theta}$ at every node and verify that it stays under a\nsmall tolerance \u2014 the trace error grows mildly with mesh size but\nnever drifts off the hoop direction.\n\n## References\n\n* Timoshenko, S. P. and Goodier, J. N. (1970) *Theory of\n Elasticity*, 3rd ed., McGraw-Hill, \u00a71.7 \u2014 principal axes,\n \u00a733 \u2014 Lam\u00e9 thick cylinder.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n *Concepts and Applications of Finite Element Analysis*, 4th ed.,\n Wiley, \u00a73.6 \u2014 stress invariants and principal stresses.\n* Sokolnikoff, I. S. (1956) *Mathematical Theory of Elasticity*,\n 2nd ed., McGraw-Hill, \u00a722 \u2014 principal-direction trajectories.\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": [ "## Build the Lam\u00e9 quarter-annulus (HEX8 EAS, plane strain)\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 = 1.0e7 # internal pressure [Pa]\nE = 2.0e11\nNU = 0.30\nN_THETA, N_RAD = 24, 16\n\ntheta = np.linspace(0.0, 0.5 * math.pi, N_THETA + 1)\nr = np.linspace(a, b, N_RAD + 1)\n\npts_list: list[list[float]] = []\nfor 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])\npts_arr = np.array(pts_list, dtype=np.float64)\n\nnx_plane = (N_THETA + 1) * (N_RAD + 1)\nn_cells = N_THETA * N_RAD\ncells = np.empty((n_cells, 9), dtype=np.int64)\ncells[:, 0] = 8\nc = 0\nfor 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\ngrid = pv.UnstructuredGrid(cells.ravel(), np.full(n_cells, 12, dtype=np.uint8), pts_arr)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": 7850.0},\n)\n\nfor 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# Internal pressure \u2014 lump as nodal forces on the bore face.\nfx_acc: dict[int, float] = {}\nfy_acc: dict[int, float] = {}\nfor 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]\nfor 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\nprint(\"Lam\u00e9 thick cylinder \u2014 principal stress + principal direction recovery\")\nprint(f\" a = {a} m, b = {b} m, p_i = {p_i / 1e6:.1f} MPa, mesh ({N_THETA}, {N_RAD})\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solve and recover the stress field\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res = m.solve_static()\nu = np.asarray(res.displacement, dtype=np.float64).ravel()\nsigma = compute_nodal_stress(m, u)\ninv = stress_invariants(sigma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Recover the principal directions per node\n\n``stress_invariants`` returns the eigenvalues only; for principal\ndirections we re-build the (n_points, 3, 3) Cauchy tensors and call\n``np.linalg.eigh`` to get sorted eigenvalues *and* eigenvectors.\n``eigh`` returns ascending eigenvalues, so column 2 is $\\sigma_{1}$.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_pts = sigma.shape[0]\nT = np.empty((n_pts, 3, 3), dtype=np.float64)\nT[:, 0, 0] = sigma[:, 0]\nT[:, 1, 1] = sigma[:, 1]\nT[:, 2, 2] = sigma[:, 2]\nT[:, 0, 1] = T[:, 1, 0] = sigma[:, 3]\nT[:, 1, 2] = T[:, 2, 1] = sigma[:, 4]\nT[:, 0, 2] = T[:, 2, 0] = sigma[:, 5]\neigvals, eigvecs = np.linalg.eigh(T)\n# eigvecs[:, :, k] is the k-th eigenvector (ascending); k=2 \u2192 \u03c31.\nv1 = eigvecs[:, :, 2]\nv3 = eigvecs[:, :, 0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verify: \u03c31 axis aligns with the analytical hoop direction\n\nAt a point $(r\\cos\\theta, r\\sin\\theta, z)$ on the annulus, the\nhoop unit vector is $\\hat{\\theta} = (-\\sin\\theta,\n\\cos\\theta, 0)$. The recovered $\\sigma_{1}$ eigenvector\nshould be parallel (or anti-parallel) to it \u2014 alignment is the\nabsolute dot product.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "r_xy = np.linalg.norm(pts_arr[:, :2], axis=1)\ntheta_node = np.arctan2(pts_arr[:, 1], pts_arr[:, 0])\nhoop_hat = np.column_stack([-np.sin(theta_node), np.cos(theta_node), np.zeros_like(theta_node)])\nalign = np.abs((v1 * hoop_hat).sum(axis=1))\n# Skip the on-axis ridge (r \u2248 0) where the hoop direction is undefined.\nmask = r_xy > 1e-6\nmean_align = float(align[mask].mean())\nworst_align = float(align[mask].min())\nprint(f\"\\n \u03c31 / hoop alignment mean |v1\u00b7\u03b8\u0302| = {mean_align:.6f} min = {worst_align:.6f}\")\nprint(f\" worst-case angle off hoop: {np.rad2deg(np.arccos(worst_align)):.3f}\u00b0\")\n\nassert mean_align > 0.999, f\"\u03c31 axis drifted from hoop direction (mean alignment {mean_align:.6f})\"\nassert worst_align > 0.95, f\"worst-case \u03c31 / hoop alignment {worst_align:.6f}\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tabulate principal stresses at the bore + outer radii\n\nClosed form for the Lam\u00e9 cylinder at any radius $r$:\n\n\\begin{align}\\sigma_\\theta(r) = p_i \\cdot \\frac{a^2\\,(b^2 + r^2)}{r^2\\,(b^2 - a^2)},\n \\qquad\n \\sigma_r(r) = -\\,p_i \\cdot \\frac{a^2\\,(b^2 - r^2)}{r^2\\,(b^2 - a^2)}.\\end{align}\n\nSo $\\sigma_1 \\approx \\sigma_\\theta$ (positive) and\n$\\sigma_3 \\approx \\sigma_r$ (negative) \u2014 at $\\theta = 0$\n/ $\\theta = \\pi/2$ the radial and hoop directions trade with the\nglobal $x / y$ axes.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print()\nprint(\n f\" {'r [m]':>7} {'\u03c31 FE [MPa]':>13} {'\u03c3_\u03b8 pub [MPa]':>14} \"\n f\"{'\u03c33 FE [MPa]':>13} {'\u03c3_r pub [MPa]':>14}\"\n)\nprint(\" \" + \"-\" * 64)\ntarget_radii = (a, 0.5 * (a + b), b)\nfor r_target in target_radii:\n # First node closest to (r_target, 0, 0).\n diff = np.linalg.norm(pts_arr[:, :2] - np.array([r_target, 0.0]), axis=1)\n i_star = int(np.argmin(diff))\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} \"\n f\"{inv['s1'][i_star] / 1e6:>11.3f} {sigma_theta_pub / 1e6:>13.3f} \"\n f\"{inv['s3'][i_star] / 1e6:>11.3f} {sigma_r_pub / 1e6:>13.3f}\"\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the \u03c31 field and principal-direction trajectories\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid_render = grid.copy()\ngrid_render.point_data[\"\u03c31 [Pa]\"] = inv[\"s1\"]\ngrid_render.point_data[\"\u03c33 [Pa]\"] = inv[\"s3\"]\ngrid_render.point_data[\"\u03c3_VM [Pa]\"] = inv[\"von_mises\"]\ngrid_render.point_data[\"v1\"] = v1\ngrid_render.point_data[\"v3\"] = v3\n\n# Sample every Mth node to keep the arrow plot readable.\nsample_stride = max(1, n_pts // 64)\nsample_idx = np.arange(0, n_pts, sample_stride)\narrow_pts = pts_arr[sample_idx]\narrow_v1 = v1[sample_idx]\n# Scale arrows uniformly to ~6 % of the bounding box.\narrow_len = 0.06 * (b - a)\n\narrow_grid = pv.PolyData(arrow_pts)\narrow_grid[\"v1\"] = arrow_v1\narrow_grid.set_active_vectors(\"v1\")\narrows = arrow_grid.glyph(orient=\"v1\", scale=False, factor=arrow_len)\n\nplotter = pv.Plotter(off_screen=True, window_size=(820, 540))\nplotter.add_mesh(\n grid_render,\n scalars=\"\u03c31 [Pa]\",\n cmap=\"plasma\",\n show_edges=False,\n scalar_bar_args={\"title\": \"\u03c31 [Pa] (max principal)\"},\n)\nplotter.add_mesh(arrows, color=\"white\", line_width=2)\nplotter.view_xy()\nplotter.camera.zoom(1.05)\nplotter.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Take-aways\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print()\nprint(\"Take-aways:\")\nprint(\n \" \u2022 stress_invariants(sigma) returns the principal *stresses* (s1 \u2265 s2 \u2265 s3) \"\n \"along with von-Mises / hydrostatic / max-shear scalars. Use it for any \"\n \"design check that needs \u03c31 (max tension) or \u03c33 (max compression).\"\n)\nprint(\n \" \u2022 Principal *directions* require the full eigendecomposition: build the \"\n \"per-node 3x3 Cauchy tensor and call np.linalg.eigh \u2014 eigvecs[:,:,2] is \"\n \"the \u03c31 axis.\"\n)\nprint(\n \" \u2022 Plotting the \u03c31 eigenvector field on top of a stress contour gives the \"\n \"stress-trajectory picture: where the load actually flows. At bore-edge \"\n \"in a Lam\u00e9 cylinder the trajectories are perfect circles.\"\n)\nprint(\n \" \u2022 Always sanity-check principal directions against an analytical reference \"\n \"where one exists (radial / hoop for cylinders, bending fibres for beams) \u2014 \"\n \"drift off the expected axis is an early warning of mesh or BC error.\"\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 }