"""Plate with a circular hole under uniaxial tension — Kirsch.
Kirsch's 1898 analytical solution for an infinite plate under
uniform remote tension :math:`\\sigma_\\infty` containing a
circular hole of radius :math:`a` gives
.. math::
\\sigma_{yy}(r, \\theta) = \\frac{\\sigma_\\infty}{2}
\\Bigl(1 + \\frac{a^{2}}{r^{2}}\\Bigr)
+ \\frac{\\sigma_\\infty}{2}
\\Bigl(1 + 3 \\frac{a^{4}}{r^{4}}\\Bigr)
\\cos(2\\theta).
At the hole boundary :math:`r = a`, the hoop stress reaches a
maximum at :math:`\\theta = \\pi/2` (perpendicular to the
remote tension axis):
.. math::
\\sigma_{yy}^{\\text{max}} = 3 \\,\\sigma_\\infty.
This is the canonical stress-concentration benchmark — the
factor of 3 is independent of ``E``, ``ν``, and the ratio
:math:`a / W` (where ``W`` is the plate width) *for an
infinite plate*. Finite plates with :math:`a / W \\lesssim 0.1`
recover :math:`K_t \\approx 3` to engineering precision.
Geometry (quarter-symmetry model)
---------------------------------
* ``W`` — plate half-width (plane stress, thickness ``t``).
* ``a`` — hole radius (:math:`a \\ll W` — we use ``W = 10 a``).
* Model the first quadrant ``x ∈ [0, W], y ∈ [0, W]`` with the
hole centred at the origin.
* Symmetry BCs: ``UX = 0`` on y-axis (``x = 0``),
``UY = 0`` on x-axis (``y = 0``).
* Load: uniform tensile traction :math:`\\sigma_\\infty` on the
far edge (``x = W``) acting in +x. Top edge (``y = W``) and
hole edge are traction-free.
References
----------
* Kirsch, E. *Die Theorie der Elastizität und die Bedürfnisse
der Festigkeitslehre.* Zeitschrift des Vereins Deutscher
Ingenieure 42 (1898), pp. 797–807.
* Timoshenko, S. P. and Goodier, J. N. *Theory of Elasticity*,
3rd ed., McGraw-Hill, 1970, §35 — standard derivation.
* Peterson, R. E. *Peterson's Stress Concentration Factors*,
3rd ed., Wiley, 2008, §4.3 — finite-width ``K_t`` correction
table (agrees with :math:`K_t = 3` at :math:`a/W \\le 0.1`).
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.3.6** plate-with-hole plane-stress family.
* **NAFEMS NL2** plate with a circular hole (linear-elastic
limit; the NAFEMS benchmark extends to plasticity).
Implementation notes
--------------------
Mapped mesh: a graded parametric grid that puts more radial
sampling near the hole (where the stress gradient is steep) and
coarser sampling out at the remote edge. Uses the
``QUAD4_PLANE`` plane-stress kernel.
The remote tension is applied as nodal forces on the far edge
(``x = W``) with trapezoid arc-length weighting — total force
= :math:`\\sigma_\\infty \\cdot W \\cdot t` on the quarter-
symmetry edge.
σ_yy at the hole edge (``x = 0, y = a``) is recovered via
forward-difference strain estimates — the framework's
``compute_nodal_stress`` helper doesn't yet cover plane-element
kernels (VERIFY-BLOCKED task #177).
"""
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 PlateWithHole(BenchmarkProblem):
"""Plate with a circular hole under uniaxial tension — Kirsch K_t = 3."""
name: str = "plate_with_hole"
description: str = (
"Quarter-symmetry plane-stress plate with a circular hole, "
"remote uniaxial tension σ_∞. Kirsch (1898) stress "
"concentration K_t = 3 at the hole edge, perpendicular "
"to the remote tension axis."
)
analysis: str = "static"
#: Hole radius [m].
a: float = 0.1
#: Plate half-width [m] (``W >> a`` for the infinite-plate limit).
W: float = 1.0
#: Out-of-plane thickness [m] (plane stress).
thickness: float = 0.01
#: Young's modulus [Pa].
E: float = 2.10e11
nu: float = 0.3
rho: float = 7850.0
#: Remote uniaxial tension [Pa].
sigma_infty: float = 1.0e7
@property
def published_values(self) -> tuple[PublishedValue, ...]:
# K_t = 3 for an infinite plate. With W = 10a the finite-
# width correction (Howland 1930) is ~0.4 %; Peterson's
# tables give K_t(a/W = 0.1) ≈ 3.04. Stay with K_t = 3
# as the Kirsch-analytical reference and use 5 % tolerance
# to soak up the finite-width correction + mesh error.
return (
PublishedValue(
name="kt_at_hole",
value=3.0 * self.sigma_infty,
unit="Pa",
source="Kirsch 1898",
formula="σ_yy_max = K_t σ_∞ with K_t = 3",
tolerance=0.10,
),
)
[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", 12))
# Parametric grid: θ ∈ [0, π/2] along the hole edge,
# s ∈ [0, 1] radially from the hole to the square outer
# boundary. Use a geometric stretch so radial nodes
# cluster near the hole where the stress gradient is
# steep.
thetas = np.linspace(0.0, np.pi / 2.0, n_theta + 1)
# Geometric spacing: s_j = (1 - q^j)/(1 - q^{n_rad}).
q = 1.25
j_idx = np.arange(n_rad + 1)
ss = (1.0 - q**j_idx) / (1.0 - q**n_rad)
nx_grid = (n_theta + 1) * (n_rad + 1)
pts = np.zeros((nx_grid, 3), dtype=np.float64)
for i, theta in enumerate(thetas):
# Inner boundary point (hole edge): (a cos θ, a sin θ).
p_in = np.array([self.a * np.cos(theta), self.a * np.sin(theta), 0.0])
# Outer boundary point along the square: intersect the
# ray (cos θ, sin θ) with the square ``max(x, y) = W``.
c = np.cos(theta)
s = np.sin(theta)
if c >= s:
# Ray exits through the right edge x = W.
t_out = self.W / c if c > 1e-12 else self.W
else:
# Ray exits through the top edge y = W.
t_out = self.W / s
p_out = np.array([t_out * c, t_out * s, 0.0])
for j, s_val in enumerate(ss):
idx = i * (n_rad + 1) + j
pts[idx] = p_in + s_val * (p_out - p_in)
# QUAD4 cells — CCW in physical xy: same pattern as NAFEMS LE1.
n_cells = n_theta * n_rad
conn = np.empty((n_cells, 5), dtype=np.int64)
conn[:, 0] = 4
c = 0
for i in range(n_theta):
for j in range(n_rad):
n00 = i * (n_rad + 1) + j
n01 = i * (n_rad + 1) + (j + 1)
n11 = (i + 1) * (n_rad + 1) + (j + 1)
n10 = (i + 1) * (n_rad + 1) + j
conn[c, 1:] = (n00, n01, n11, n10)
c += 1
grid = pv.UnstructuredGrid(
conn.ravel(),
np.full(n_cells, 9, dtype=np.int64), # VTK_QUAD
pts,
)
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.QUAD4_PLANE,
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
real=(self.thickness,),
)
# Symmetry BCs.
for k, p in enumerate(pts):
if p[1] < 1e-9: # x-axis
m.fix(nodes=int(k + 1), dof="UY", value=0.0)
if p[0] < 1e-9: # y-axis
m.fix(nodes=int(k + 1), dof="UX", value=0.0)
# Remote tension on the far edge ``x = W`` — apply +x
# nodal forces with trapezoid arc-length weighting. The
# right edge is the subset of nodes at j = n_rad with
# θ in the range where the ray exits through x = W
# (``θ ∈ [0, π/4]``). Same for the top edge y = W at
# ``θ ∈ [π/4, π/2]``. We apply tension on the right edge
# only — symmetric tension on both x-facing edges doubles
# the far-field σ_∞ counting.
#
# For a quarter-symmetry plate under uniaxial σ_∞ in the
# x-direction, the boundary traction is σ_∞ · n̂ = σ_∞ x̂
# on the right edge only; the top edge is traction-free.
#
# Total applied force on the right edge = σ_∞ · W · t
# (per-quadrant-edge tension = σ_∞ × plate half-width × thickness).
right_nodes: list[int] = []
for k, p in enumerate(pts):
if abs(p[0] - self.W) < 1e-9:
right_nodes.append(k)
# Sort by y (so adjacent nodes form edge segments).
right_nodes.sort(key=lambda k: pts[k, 1])
fx_by_node: dict[int, float] = {}
for seg in range(len(right_nodes) - 1):
a_idx = right_nodes[seg]
b_idx = right_nodes[seg + 1]
ds = float(np.abs(pts[b_idx, 1] - pts[a_idx, 1]))
F_seg = self.sigma_infty * ds * self.thickness
fx_by_node[a_idx] = fx_by_node.get(a_idx, 0.0) + 0.5 * F_seg
fx_by_node[b_idx] = fx_by_node.get(b_idx, 0.0) + 0.5 * F_seg
for n_idx, fx in fx_by_node.items():
m.apply_force(int(n_idx + 1), fx=fx)
# Stash the hole-top node (θ=π/2, j=0) for extract. This
# is the point where σ_yy peaks.
hole_top_idx = n_theta * (n_rad + 1) + 0
assert np.isclose(pts[hole_top_idx, 0], 0.0)
assert np.isclose(pts[hole_top_idx, 1], self.a)
m._bench_hole_top_idx = hole_top_idx # type: ignore[attr-defined]
m._bench_pts = pts # type: ignore[attr-defined]
m._bench_n_rad = n_rad # type: ignore[attr-defined]
return m