"""Thick-walled cylinder under internal pressure — Lamé problem.
Classical pressure-vessel benchmark. A long cylinder with inner
radius :math:`a`, outer radius :math:`b`, subjected to internal
pressure :math:`p_i` (no external pressure) deforms with pure
radial displacement ``u_r(r)`` and develops radial / hoop
stresses that Lamé derived in 1852.
Plane-strain approximation: the cylinder is axially constrained
(``\\varepsilon_z = 0``). We discretise a thin 3D slab with
UZ pinned on both faces, so the solid-element assembly produces
the plane-strain stiffness without needing a dedicated
axisymmetric kernel.
Geometry: quarter-annulus model (``x ≥ 0, y ≥ 0``) with
symmetry BCs on the cut planes.
Closed-form quantities
----------------------
At radius :math:`r \\in [a, b]`:
.. math::
\\sigma_r(r) &= \\frac{p_i a^{2}}{b^{2}-a^{2}}
\\Bigl(1 - \\frac{b^{2}}{r^{2}}\\Bigr), \\\\
\\sigma_\\theta(r) &= \\frac{p_i a^{2}}{b^{2}-a^{2}}
\\Bigl(1 + \\frac{b^{2}}{r^{2}}\\Bigr), \\\\
u_r(r) &= \\frac{p_i\\,a^{2}\\,r}{E(b^{2}-a^{2})}
\\Bigl[(1 - \\nu - 2\\nu^{2})
+ \\frac{b^{2}(1+\\nu)}{r^{2}}\\Bigr]
\\quad \\text{(plane strain).}
Notable published values:
* :math:`\\sigma_\\theta(a) = p_i (a^{2} + b^{2})/(b^{2}-a^{2})`
— hoop stress concentration at the inner surface.
* :math:`\\sigma_r(a) = -p_i` — radial stress at the inner
surface equals the applied pressure (sign convention:
compression).
References
----------
* Lamé, G. *Leçons sur la Théorie Mathématique de l'Élasticité
des Corps Solides*, Bachelier, 1852 — original derivation.
* Timoshenko, S. P. and Goodier, J. N. *Theory of Elasticity*,
3rd ed., McGraw-Hill, 1970, §33 — standard textbook
exposition of the Lamé problem.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.1.14** thick-walled cylinder family.
* **NAFEMS FV-?** (cylinder-pressure benchmarks in the NAFEMS
pressure-vessel tests).
Implementation notes
--------------------
Mapped quarter-annulus mesh: regular ``(θ, s)`` grid with
``θ ∈ [0, π/2]`` and ``s ∈ [0, 1]`` interpolating radially from
``r = a`` to ``r = b``. One axial cell layer with UZ = 0 on
both faces (plane strain).
Internal pressure lumped onto the inner-surface nodes via
trapezoid arc-length weighting; nodal accumulation is done
locally before calling ``Model.apply_force`` to avoid the
assignment-not-accumulation semantics of ``_f_impl``.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any
import numpy as np
import pyvista as pv
import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.validation._benchmark import (
BenchmarkProblem,
PublishedValue,
)
[docs]
@dataclass
class LameCylinder(BenchmarkProblem):
"""Thick-walled cylinder under internal pressure (plane strain)."""
name: str = "lame_cylinder"
description: str = (
"Thick-walled cylinder under internal pressure — Lamé "
"closed-form radial displacement + hoop stress at r = a "
"(Timoshenko & Goodier 1970 §33)."
)
analysis: str = "static"
#: Inner radius [m].
a: float = 0.01
#: Outer radius [m].
b: float = 0.02
#: Axial thickness of the 3D slab [m] (plane-strain model).
t_axial: float = 0.001
#: Young's modulus [Pa].
E: float = 2.10e11
nu: float = 0.3
rho: float = 7850.0
#: Internal pressure [Pa].
p_i: float = 1.0e8
@property
def published_values(self) -> tuple[PublishedValue, ...]:
a2 = self.a * self.a
b2 = self.b * self.b
# Hoop stress at the inner surface.
sigma_theta_a = self.p_i * (a2 + b2) / (b2 - a2)
# Radial displacement at the inner surface (plane strain).
coef = self.p_i * a2 * self.a / (self.E * (b2 - a2))
ur_a = coef * ((1 - self.nu - 2 * self.nu**2) + b2 * (1 + self.nu) / a2)
# Hoop stress at the outer surface.
sigma_theta_b = self.p_i * (2 * a2) / (b2 - a2)
return (
PublishedValue(
name="ur_at_a",
value=ur_a,
unit="m",
source="Timoshenko & Goodier 1970 §33",
formula=("u_r(a) = p_i a^2 a / (E (b^2 - a^2)) * ((1-ν-2ν^2) + b^2(1+ν)/a^2)"),
tolerance=0.03,
),
PublishedValue(
name="sigma_theta_at_a",
value=sigma_theta_a,
unit="Pa",
source="Timoshenko & Goodier 1970 §33",
formula="σ_θ(a) = p_i (a^2 + b^2) / (b^2 - a^2)",
tolerance=0.08,
),
PublishedValue(
name="sigma_theta_at_b",
value=sigma_theta_b,
unit="Pa",
source="Timoshenko & Goodier 1970 §33",
formula="σ_θ(b) = 2 p_i a^2 / (b^2 - a^2)",
tolerance=0.08,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
n_theta = int(mesh_params.get("n_theta", 16))
n_rad = int(mesh_params.get("n_rad", 6))
thetas = np.linspace(0.0, np.pi / 2.0, n_theta + 1)
ss = np.linspace(0.0, 1.0, n_rad + 1)
# Two z-layers (top + bottom) for the plane-strain slab.
zs = np.array([0.0, self.t_axial])
nx_plane = (n_theta + 1) * (n_rad + 1)
pts = np.zeros((nx_plane * 2, 3), dtype=np.float64)
for kz, z in enumerate(zs):
for i, theta in enumerate(thetas):
for j, s in enumerate(ss):
r = self.a + s * (self.b - self.a)
idx = kz * nx_plane + i * (n_rad + 1) + j
pts[idx] = (r * np.cos(theta), r * np.sin(theta), z)
# Hex8 cells (one axial layer × n_theta × n_rad angular/radial).
n_cells = n_theta * n_rad
cell_connectivity = np.empty((n_cells, 9), dtype=np.int64)
cell_connectivity[:, 0] = 8
c = 0
for i in range(n_theta):
for j in range(n_rad):
# Bottom-face (z=0) corners — CCW in xy.
n00_b = 0 * nx_plane + i * (n_rad + 1) + j
n10_b = 0 * nx_plane + i * (n_rad + 1) + (j + 1)
n11_b = 0 * nx_plane + (i + 1) * (n_rad + 1) + (j + 1)
n01_b = 0 * nx_plane + (i + 1) * (n_rad + 1) + j
# Top-face corners — same ordering, offset in z.
n00_t = n00_b + nx_plane
n10_t = n10_b + nx_plane
n11_t = n11_b + nx_plane
n01_t = n01_b + nx_plane
# VTK_HEXAHEDRON ordering: bottom CCW, top CCW.
cell_connectivity[c, 1:] = (
n00_b,
n10_b,
n11_b,
n01_b,
n00_t,
n10_t,
n11_t,
n01_t,
)
c += 1
grid = pv.UnstructuredGrid(
cell_connectivity.ravel(),
np.full(n_cells, 12, dtype=np.int64), # VTK_HEXAHEDRON
pts,
)
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.HEX8,
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
)
# Symmetry + plane-strain BCs.
for k, p in enumerate(pts):
if p[1] < 1e-9: # y = 0 (x-axis symmetry plane)
m.fix(nodes=int(k + 1), dof="UY", value=0.0)
if p[0] < 1e-9: # x = 0 (y-axis symmetry plane)
m.fix(nodes=int(k + 1), dof="UX", value=0.0)
# Plane strain: pin UZ on every node.
m.fix(nodes=int(k + 1), dof="UZ", value=0.0)
# Internal pressure on r = a. Inner nodes are at (θ_i,
# s=0) on both z layers. Accumulate contributions locally.
fx_by_node: dict[int, float] = {}
fy_by_node: dict[int, float] = {}
for kz in range(2):
inner_idx = [kz * nx_plane + i * (n_rad + 1) + 0 for i in range(n_theta + 1)]
for segment in range(n_theta):
a_idx = inner_idx[segment]
b_idx = inner_idx[segment + 1]
pa_pt = pts[a_idx]
pb_pt = pts[b_idx]
ds = float(np.linalg.norm(pb_pt - pa_pt))
# Outward normal from the cavity surface into the
# material — points radially outward from the cyl.
# axis (away from origin in xy).
mid = 0.5 * (pa_pt + pb_pt)
r_xy = np.array([mid[0], mid[1], 0.0])
nrm = np.linalg.norm(r_xy)
if nrm < 1e-12:
continue
outward = r_xy / nrm
# Force per unit axial length = pressure × arc
# length. Each z-layer covers half the axial
# thickness (we distribute the slab load as
# p × ds × (t_axial / 2) to each z-layer's node
# pair — trapezoid rule in z).
F_seg = self.p_i * ds * (self.t_axial / 2.0)
for n_idx in (a_idx, b_idx):
fx_by_node[n_idx] = fx_by_node.get(n_idx, 0.0) + 0.5 * F_seg * outward[0]
fy_by_node[n_idx] = fy_by_node.get(n_idx, 0.0) + 0.5 * F_seg * outward[1]
for n_idx, fx in fx_by_node.items():
fy = fy_by_node[n_idx]
m.apply_force(int(n_idx + 1), fx=fx, fy=fy)
# Stash inner / outer-surface node id rows (x-axis, y=0,
# z=0) for extract().
inner_x_axis = 0 * nx_plane + 0 * (n_rad + 1) + 0
outer_x_axis = 0 * nx_plane + 0 * (n_rad + 1) + n_rad
assert np.isclose(pts[inner_x_axis, 0], self.a)
assert np.isclose(pts[outer_x_axis, 0], self.b)
m._bench_inner_x_idx = inner_x_axis # type: ignore[attr-defined]
m._bench_outer_x_idx = outer_x_axis # type: ignore[attr-defined]
m._bench_pts = pts # type: ignore[attr-defined]
m._bench_n_theta = n_theta # type: ignore[attr-defined]
m._bench_n_rad = n_rad # type: ignore[attr-defined]
m._bench_nx_plane = nx_plane # type: ignore[attr-defined]
return m