.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_lame_cylinder.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_verification_example_verify_lame_cylinder.py: Lamé thick-walled cylinder under internal pressure ================================================== Classical pressure-vessel benchmark: a long cylinder with inner radius :math:`a`, outer radius :math:`b`, subjected to internal pressure :math:`p_i` and held at zero axial strain (plane strain). The closed-form Lamé (1852) solution gives purely radial displacement :math:`u_r(r)` and the radial / hoop stresses .. math:: \sigma_r(r) = \frac{p_i a^{2}}{b^{2} - a^{2}} \left(1 - \frac{b^{2}}{r^{2}}\right), .. math:: \sigma_\theta(r) = \frac{p_i a^{2}}{b^{2} - a^{2}} \left(1 + \frac{b^{2}}{r^{2}}\right), .. math:: u_r(r) = \frac{p_i a^{2} r}{E (b^{2} - a^{2})} \left[(1 - \nu - 2 \nu^{2}) + \frac{b^{2} (1 + \nu)}{r^{2}}\right] \quad (\text{plane strain}). Two textbook quantities make the cleanest comparison: * :math:`u_r(a)` — radial displacement at the bore. * :math:`\sigma_\theta(a) = p_i (a^{2} + b^{2}) / (b^{2} - a^{2})` — hoop stress concentration at the bore (always greater than :math:`p_i` for any thick wall :math:`b > a`). Implementation -------------- Quarter-annulus model (:math:`x \ge 0,\ y \ge 0`) meshed as a single axial slab of HEX8 elements with ``KEYOPT(1)=1`` enhanced strain to suppress shear locking near the bore. Symmetry / plane-strain BCs: * :math:`u_x = 0` on the cut plane :math:`x = 0`. * :math:`u_y = 0` on the cut plane :math:`y = 0`. * :math:`u_z = 0` on every node (plane strain). Internal pressure is applied as a consistent set of nodal forces on the inner surface (one trapezoid per circumferential segment, half-load per axial layer). The hoop stress is read out on the bore from the small-strain finite-difference of the solved displacement field (only the average across the bore is needed for the partition-of-unity-style verification check). References ---------- * Lamé, G. (1852) *Leçons sur la Théorie Mathématique de l'Élasticité des Corps Solides*, Bachelier — original derivation. * Timoshenko, S. P. and Goodier, J. N. (1970) *Theory of Elasticity*, 3rd ed., McGraw-Hill, §33. * Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for Stress and Strain*, 6th ed., McGraw-Hill, Table 32 case 1b (thick-wall cylinder, internal pressure). Cross-references (problem-id only; no vendor content vendored): * (Reference benchmark — the canonical Lamé thick-cylinder citation): Lamé, G., 1852 — *Leçons sur la théorie de l'élasticité*: stresses in a long cylinder under internal pressure*. .. GENERATED FROM PYTHON SOURCE LINES 73-84 .. code-block:: Python from __future__ import annotations import math import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS .. GENERATED FROM PYTHON SOURCE LINES 85-90 Problem data ------------ Steel-like elastic constants; pressure modest enough that the linear-elastic regime is appropriate. .. GENERATED FROM PYTHON SOURCE LINES 90-111 .. code-block:: Python a = 0.10 # bore radius [m] b = 0.20 # outer radius [m] t_axial = 0.02 # axial slab thickness [m] (purely geometric — plane strain) p_i = 1.0e7 # internal pressure [Pa] E = 2.0e11 # Young's modulus [Pa] NU = 0.30 # Poisson's ratio RHO = 7850.0 # density [kg/m^3] # Closed-form quantities ------------------------------------------------- ur_a_published = ( p_i * a**2 * a / (E * (b**2 - a**2)) * ((1 - NU - 2 * NU**2) + b**2 * (1 + NU) / a**2) ) sigma_theta_a_published = p_i * (a**2 + b**2) / (b**2 - a**2) print(f"problem: a={a} m, b={b} m, p_i={p_i:.1e} Pa, E={E:.1e} Pa, nu={NU}") print(f"u_r(a) = {ur_a_published:.6e} m (Lamé)") print(f"sigma_theta(a) = {sigma_theta_a_published:.6e} Pa") print(f"sigma_theta(a) / p_i = {sigma_theta_a_published / p_i:.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none problem: a=0.1 m, b=0.2 m, p_i=1.0e+07 Pa, E=2.0e+11 Pa, nu=0.3 u_r(a) = 9.533333e-06 m (Lamé) sigma_theta(a) = 1.666667e+07 Pa sigma_theta(a) / p_i = 1.6667 .. GENERATED FROM PYTHON SOURCE LINES 112-120 Build a quarter-annulus HEX8 mesh --------------------------------- A single axial slab (one layer in :math:`z`) carries the plane-strain solution since UZ is pinned everywhere. The circumferential / radial discretisation is moderate so the docs build runs fast; convergence vs published values is studied in the validation suite (:class:`femorph_solver.validation.problems.LameCylinder`). .. GENERATED FROM PYTHON SOURCE LINES 120-162 .. code-block:: Python n_theta = 12 n_rad = 8 theta = np.linspace(0.0, 0.5 * math.pi, n_theta + 1) r = np.linspace(a, b, n_rad + 1) pts: list[list[float]] = [] for kz in (0.0, t_axial): for ti in theta: for rj in r: pts.append([rj * math.cos(ti), rj * math.sin(ti), kz]) pts_arr = np.array(pts, dtype=np.float64) nx_plane = (n_theta + 1) * (n_rad + 1) n_cells = n_theta * n_rad cells = np.empty((n_cells, 9), dtype=np.int64) cells[:, 0] = 8 c = 0 for i in range(n_theta): for j in range(n_rad): n00b = i * (n_rad + 1) + j n10b = i * (n_rad + 1) + (j + 1) n11b = (i + 1) * (n_rad + 1) + (j + 1) n01b = (i + 1) * (n_rad + 1) + j cells[c, 1:] = ( n00b, n10b, n11b, n01b, n00b + nx_plane, n10b + nx_plane, n11b + nx_plane, n01b + nx_plane, ) c += 1 grid = pv.UnstructuredGrid( cells.ravel(), np.full(n_cells, 12, dtype=np.uint8), # VTK_HEXAHEDRON pts_arr, ) print(f"mesh: {grid.n_points} nodes, {grid.n_cells} HEX8 cells") .. rst-class:: sphx-glr-script-out .. code-block:: none mesh: 234 nodes, 96 HEX8 cells .. GENERATED FROM PYTHON SOURCE LINES 163-165 Wrap in a femorph-solver Model and stamp BCs -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 165-180 .. code-block:: Python m = femorph_solver.Model.from_grid(grid) m.assign( ELEMENTS.HEX8(integration="enhanced_strain"), material={"EX": E, "PRXY": NU, "DENS": RHO}, ) # Symmetry + plane-strain BCs. for k, p in enumerate(pts_arr): if p[0] < 1e-9: m.fix(nodes=[int(k + 1)], dof="UX") if p[1] < 1e-9: m.fix(nodes=[int(k + 1)], dof="UY") m.fix(nodes=[int(k + 1)], dof="UZ") .. GENERATED FROM PYTHON SOURCE LINES 181-189 Apply internal pressure as a consistent nodal-force ring -------------------------------------------------------- Each :math:`(r = a)` segment between two circumferential nodes carries a force :math:`F = p_i \cdot \mathrm{ds} \cdot (t/2)` on each axial layer (trapezoid rule in :math:`z`); the half-load splits between the two endpoints; the direction is the outward normal at the segment midpoint. .. GENERATED FROM PYTHON SOURCE LINES 189-209 .. code-block:: Python fx_acc: dict[int, float] = {} fy_acc: dict[int, float] = {} for kz in (0, 1): inner = [kz * nx_plane + i * (n_rad + 1) + 0 for i in range(n_theta + 1)] for seg in range(n_theta): ai, bi = inner[seg], inner[seg + 1] ds = float(np.linalg.norm(pts_arr[bi] - pts_arr[ai])) mid = 0.5 * (pts_arr[ai] + pts_arr[bi]) rxy = np.array([mid[0], mid[1]]) nrm = float(np.linalg.norm(rxy)) outward = rxy / nrm if nrm > 1e-12 else np.zeros(2) F_seg = p_i * ds * (t_axial / 2.0) for n_idx in (ai, bi): fx_acc[n_idx] = fx_acc.get(n_idx, 0.0) + 0.5 * F_seg * outward[0] fy_acc[n_idx] = fy_acc.get(n_idx, 0.0) + 0.5 * F_seg * outward[1] for n_idx, fx in fx_acc.items(): fy = fy_acc[n_idx] m.apply_force(int(n_idx + 1), fx=fx, fy=fy) .. GENERATED FROM PYTHON SOURCE LINES 210-212 Solve and recover the bore displacement --------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 212-231 .. code-block:: Python res = m.solve() u_flat = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 3) # Inner-bore node on the +x-axis (i=0, j=0, kz=0). By symmetry # its UX equals u_r(a). inner_x = 0 ur_a_computed = float(u_flat[inner_x, 0]) err_ur = (ur_a_computed - ur_a_published) / ur_a_published print() print(f"u_r(a) computed = {ur_a_computed:.6e} m") print(f"u_r(a) published = {ur_a_published:.6e} m") print(f"relative error = {err_ur * 100:+.3f} %") # 5 % tolerance — moderate-mesh HEX8 EAS captures u_r(a) to a # couple percent on this geometry (n_theta=12, n_rad=8); the # validation suite drives convergence to < 1 % at finer meshes. assert abs(err_ur) < 0.05, f"u_r(a) error {err_ur:.2%} exceeds 5%" .. rst-class:: sphx-glr-script-out .. code-block:: none u_r(a) computed = 9.513374e-06 m u_r(a) published = 9.533333e-06 m relative error = -0.209 % .. GENERATED FROM PYTHON SOURCE LINES 232-241 Estimate sigma_theta(a) by central-difference of UY in theta ------------------------------------------------------------ At the y=0 cut plane :math:`u_\theta(r, 0) = 0`, so :math:`\varepsilon_\theta \approx u_y(r, \Delta\theta) / r\Delta\theta` at the next circumferential node. Combine with the radial strain :math:`\varepsilon_r = \partial u_x / \partial x` from the next radial node, then plug into the plane-strain constitutive law. .. GENERATED FROM PYTHON SOURCE LINES 241-262 .. code-block:: Python idx_theta_next = 0 * nx_plane + 1 * (n_rad + 1) + 0 # i=1, j=0, kz=0 idx_rad_next = 0 * nx_plane + 0 * (n_rad + 1) + 1 # i=0, j=1, kz=0 y_next = float(pts_arr[idx_theta_next, 1]) eps_theta = float(u_flat[idx_theta_next, 1]) / y_next dx = float(pts_arr[idx_rad_next, 0] - pts_arr[inner_x, 0]) eps_r = float(u_flat[idx_rad_next, 0] - u_flat[inner_x, 0]) / dx Cfac = E / ((1 + NU) * (1 - 2 * NU)) sigma_theta_a_computed = Cfac * ((1 - NU) * eps_theta + NU * eps_r) err_sigma = (sigma_theta_a_computed - sigma_theta_a_published) / sigma_theta_a_published print() print(f"sigma_theta(a) computed = {sigma_theta_a_computed:.4e} Pa") print(f"sigma_theta(a) published = {sigma_theta_a_published:.4e} Pa") print(f"relative error = {err_sigma * 100:+.3f} %") # 10 % tolerance — the finite-difference stress recovery on a # coarse HEX8 mesh is intrinsically noisier than the displacement # match; full-field stress recovery on a refined mesh closes the # gap to <2 % in the validation suite. assert abs(err_sigma) < 0.10, f"sigma_theta(a) error {err_sigma:.2%} exceeds 10%" .. rst-class:: sphx-glr-script-out .. code-block:: none sigma_theta(a) computed = 1.7750e+07 Pa sigma_theta(a) published = 1.6667e+07 Pa relative error = +6.499 % .. GENERATED FROM PYTHON SOURCE LINES 263-268 Render the radial displacement field on the deformed quarter-annulus -------------------------------------------------------------------- Scatter the DOF-indexed displacement onto :math:`(n, 3)` point data, then warp the mesh and colour by :math:`u_r`. .. GENERATED FROM PYTHON SOURCE LINES 268-297 .. code-block:: Python grid_plot = femorph_solver.io.static_result_to_grid(m, res) ux = grid_plot.point_data["displacement"][:, 0] uy = grid_plot.point_data["displacement"][:, 1] xy = np.asarray(grid_plot.points)[:, :2] r_node = np.linalg.norm(xy, axis=1) ur = (ux * xy[:, 0] + uy * xy[:, 1]) / np.maximum(r_node, 1e-12) grid_plot.point_data["u_r"] = ur warped = grid_plot.warp_by_vector("displacement", factor=2.0e3) plotter = pv.Plotter(off_screen=True, window_size=(720, 480)) plotter.add_mesh( grid_plot, style="wireframe", color="grey", opacity=0.35, label="undeformed", ) plotter.add_mesh( warped, scalars="u_r", cmap="plasma", show_edges=True, scalar_bar_args={"title": "u_r [m] (deformed ×2000)"}, ) plotter.view_xy() plotter.camera.zoom(1.1) plotter.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_lame_cylinder_001.png :alt: example verify lame cylinder :srcset: /gallery/verification/images/sphx_glr_example_verify_lame_cylinder_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.231 seconds) .. _sphx_glr_download_gallery_verification_example_verify_lame_cylinder.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_verify_lame_cylinder.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_lame_cylinder.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_lame_cylinder.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_