r""" 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). Vendor cross-references ----------------------- ====================================== ========================= ============================================ Source Reported σ_θ(a) [MPa] Problem ID / location ====================================== ========================= ============================================ Closed form (Lamé 1852) 166.67 Timoshenko & Goodier 1970 §33 Timoshenko & Goodier (1970) 166.67 Theory of Elasticity §33 MAPDL Verification Manual ≈ 166.7 VM-25 Stresses in a long cylinder Abaqus Verification Manual ≈ 166.7 AVM 1.1.14 thick-walled cylinder family ====================================== ========================= ============================================ """ from __future__ import annotations import math import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS # %% # Problem data # ------------ # # Steel-like elastic constants; pressure modest enough that the # linear-elastic regime is appropriate. 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}") # %% # 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`). 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") # %% # Wrap in a femorph-solver Model and stamp BCs # -------------------------------------------- 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") # %% # 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. 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) # %% # Solve and recover the bore displacement # --------------------------------------- res = m.solve_static() 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%" # %% # 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. 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%" # %% # 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`. 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()