Source code for femorph_solver.validation._problems.pinched_ring

"""Pinched ring — diametrical deflection under two opposed point loads.

Classical thin-ring benchmark.  A circular ring of mean radius
:math:`R`, bending stiffness :math:`EI`, loaded by two equal and
opposite point forces :math:`P` acting along a diameter.
Castigliano's theorem gives the diametrical deflection

.. math::

    \\delta_{\\text{diam}} = \\frac{P R^{3}}{E I}
        \\Bigl(\\frac{\\pi}{4} - \\frac{2}{\\pi}\\Bigr)
    \\approx 0.14868 \\, \\frac{P R^{3}}{E I},

measuring the change in ring diameter along the line of action
of the loads (the loaded points move toward each other by
:math:`\\delta_{\\text{diam}} / 2` each).

References
----------
* Timoshenko, S. P. and Young, D. H.  *Elements of Strength of
  Materials*, 5th ed., Van Nostrand, 1968, §79 — derivation of
  the pinched-ring deflection using Castigliano's theorem.
* Roark's *Formulas for Stress and Strain*, 8th ed., 2011,
  Table 9.2 Case 13 — tabulated pinched-ring formulas.

Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):

* **Abaqus AVM 1.5.x** pinched-ring family.
* **NAFEMS R0011** pinched-ring benchmark (quarter-symmetry
  variant).

Implementation notes
--------------------

Quarter-symmetry model — one quadrant of the ring
(:math:`\\theta \\in [0, \\pi/2]`) modelled with a chain of
``BEAM2`` elements.  Beam local axes require the element
tangent direction + a reference orientation; we build the
circular chord with nodes on the ring and let the BEAM2
kernel's auto-orientation pick a consistent transverse axis.

BCs: symmetry about the x-axis (``UY = ROTX = ROTZ = 0`` at
the loaded point ``(R, 0)``) and symmetry about the y-axis
(``UX = ROTY = ROTZ = 0`` at the free point ``(0, R)`` —
actually for the quarter model we fix the point where the
diameter exits, not the loaded point).

We apply half of the load (``P/2``) as a nodal force at the
loaded symmetry node.  By symmetry, the quarter-ring
deflection along the line of action equals
:math:`\\delta_{\\text{diam}} / 2`.
"""

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 PinchedRing(BenchmarkProblem): """Pinched-ring diametrical deflection (quarter-symmetry).""" name: str = "pinched_ring" description: str = ( "Thin circular ring under two opposed point loads — " "Castigliano closed-form diametrical deflection " "(Timoshenko & Young 1968 §79)." ) analysis: str = "static" #: Ring mean radius [m]. R: float = 0.1 #: Cross-section width [m]. width: float = 0.01 #: Cross-section height (radial-thickness) [m]. height: float = 0.005 #: Young's modulus [Pa]. E: float = 2.0e11 nu: float = 0.3 rho: float = 7850.0 #: Magnitude of each applied point load [N]. P: float = 10.0 @property def published_values(self) -> tuple[PublishedValue, ...]: I = self.width * self.height**3 / 12.0 # noqa: E741 # The diametrical deflection (full-ring) equals P·R³/(E·I) · (π/4 - 2/π). # In the quarter-symmetry model the loaded-point UX displacement is # half of that (each end moves by δ/2 along its own symmetry axis). delta = self.P * self.R**3 / (self.E * I) * (np.pi / 4.0 - 2.0 / np.pi) return ( PublishedValue( name="loaded_point_ux", value=delta / 2.0, unit="m", source="Timoshenko & Young 1968 §79", formula="u_x(loaded) = (P R^3 / (E I)) (π/4 - 2/π) / 2", tolerance=0.05, ), )
[docs] def build_model(self, **mesh_params: Any) -> femorph_solver.Model: n_elem = int(mesh_params.get("n_elem", 40)) # Quarter arc parameterised by ``θ ∈ [0, π/2]``. Node 0 # is at (R, 0) — the loaded point. Node n_elem is at # (0, R) — the "pinched" symmetry point on the y-axis. thetas = np.linspace(0.0, np.pi / 2.0, n_elem + 1) pts = np.column_stack( (self.R * np.cos(thetas), self.R * np.sin(thetas), np.zeros_like(thetas)) ) cells = np.empty((n_elem, 3), dtype=np.int64) cells[:, 0] = 2 cells[:, 1] = np.arange(n_elem) cells[:, 2] = np.arange(1, n_elem + 1) grid = pv.UnstructuredGrid( cells.ravel(), np.full(n_elem, 3, dtype=np.int64), pts, ) m = femorph_solver.Model.from_grid(grid) A = self.width * self.height Iz = self.width * self.height**3 / 12.0 Iy = self.height * self.width**3 / 12.0 J = (1.0 / 3.0) * min(self.width, self.height) ** 3 * max(self.width, self.height) m.assign( ELEMENTS.BEAM2, material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho}, real=(A, Iz, Iy, J), ) # Symmetry BCs on the two cut ends of the quarter arc. # Loaded point (R, 0): y-axis of symmetry is perpendicular # to the load → UY = 0 (no transverse displacement across # the load plane), rotations about the in-plane axes also # zero (ROTZ = 0 for in-plane bending symmetry; ROTX, # ROTY pinned to kill spurious 3D modes). m.fix(nodes=1, dof="UY", value=0.0) m.fix(nodes=1, dof="UZ", value=0.0) m.fix(nodes=1, dof="ROTX", value=0.0) m.fix(nodes=1, dof="ROTY", value=0.0) m.fix(nodes=1, dof="ROTZ", value=0.0) # Pinched symmetry point (0, R): UX = 0 along the # y-symmetry cut; ROTZ = 0 to enforce full-ring symmetry # (the cross-section stays aligned with the radial # direction at the cut); UZ and the out-of-plane rotations # pinned to suppress spurious 3D modes. m.fix(nodes=n_elem + 1, dof="UX", value=0.0) m.fix(nodes=n_elem + 1, dof="UZ", value=0.0) m.fix(nodes=n_elem + 1, dof="ROTX", value=0.0) m.fix(nodes=n_elem + 1, dof="ROTY", value=0.0) m.fix(nodes=n_elem + 1, dof="ROTZ", value=0.0) # Apply P/2 at the loaded point along +x (away from the # origin — the symmetric full-ring load pair is "squeeze # inward" which in the quarter-model maps to a force # component outward along the x-axis for the quarter-arc). # Actually for a pinched ring compressed along x, the # loaded point moves in -x. In the quarter-symmetry # model the applied load at the loaded point is +P/2 # inward (i.e. -x), but by Castigliano's derivation we # apply P/2 inward toward the origin to match the # "compression" sense of the classical problem. m.apply_force(1, fx=-self.P / 2.0) m._bench_loaded_node_id = 1 # type: ignore[attr-defined] return m
[docs] def extract(self, model: femorph_solver.Model, result: Any, name: str) -> float: if name != "loaded_point_ux": raise KeyError(f"unknown quantity {name!r}") u = np.asarray(result.displacement, dtype=np.float64) n_nodes = model.grid.n_points dofs_per_node = u.size // n_nodes u2 = u.reshape(n_nodes, dofs_per_node) loaded = model._bench_loaded_node_id - 1 # type: ignore[attr-defined] # Report the magnitude of the inward displacement at the # loaded point (the loaded point moves toward the origin # by δ/2, so u_x is negative for a load pushing -x; return # the absolute value to compare against the positive # published deflection). return float(-u2[loaded, 0])