Plate with a circular hole — Kirsch stress concentration

Classical 2D elasticity benchmark: a thin rectangular plate with a small circular hole, loaded by uniform remote tension \(\sigma_{\infty}\) along the \(x\)-axis. Kirsch’s 1898 closed form for the infinite plate gives a stress field that, on the hole edge \(r = a\), peaks at the points perpendicular to the loading direction (\(\theta = \pm\pi/2\)):

(1)\[K_{t} \;\equiv\; \frac{\sigma_{yy}^{(\max)}}{\sigma_{\infty}} \;=\; 3 \qquad \text{(infinite plate, hole radius } a \text{)}.\]

This is the textbook 3:1 stress concentration, the foundational example for fatigue-design notch factors and the canonical first-principles check for stress-recovery correctness.

Implementation

Quarter-symmetry plane-stress model (PlateWithHole) on a \(16 \times 12\) QUAD4 mapped mesh that radially clusters nodes near the hole where the stress gradient is steep. Symmetry pins \(u_y = 0\) on \(y = 0\) and \(u_x = 0\) on \(x = 0\); the right edge (\(x = W/2\)) carries a consistent traction load equivalent to a uniform \(\sigma_{\infty}\) remote tension. Finite-width correction (Howland 1930; Peterson Table 4-1) gives \(K_t \approx 3.04\) at \(a/W = 0.1\), so the expected FE result lies between the Kirsch infinite-plate value (3.0) and the finite-width correction (~3.04) plus a small mesh-discretisation overhead.

References

• Kirsch, G. (1898) “Die Theorie der Elastizität und die Bedürfnisse der Festigkeitslehre,” Zeitschrift des Vereins Deutscher Ingenieure 42, 797–807 — the original closed-form derivation.

• Timoshenko, S. P. and Goodier, J. N. (1970) Theory of Elasticity, 3rd ed., McGraw-Hill, §35 (plate with a circular hole).

• Howland, R. C. J. (1930) “On the stresses in the neighbourhood of a circular hole in a strip under tension,” Phil. Trans. R. Soc. A 229, 49–86 — finite- width correction.

• Peterson, R. E. (1974) Stress Concentration Factors, Wiley, Table 4-1 (uniaxial stress, hole in finite plate).

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §5.5 (stress concentration benchmark).

Cross-references (problem-id only):

• (Reference benchmark — the canonical Kirsch infinite-plate citation): Kirsch, G., 1898 — VDI 42: stress concentration around a hole in a plate.*

• Abaqus AVM 1.3.6 elliptic-membrane / plate-with-hole.

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation.problems import PlateWithHole

Build the model from the validation problem class

problem = PlateWithHole()
m = problem.build_model()

print(
    f"plate-with-hole mesh: {m.grid.n_points} nodes, {m.grid.n_cells} QUAD4 cells; "
    f"a = {problem.a} m, W = {problem.W} m, σ_∞ = {problem.sigma_infty / 1e6:.1f} MPa"
)
print(f"E = {problem.E / 1e9:.0f} GPa, ν = {problem.nu}, a/W = {problem.a / problem.W:.2f}")

plate-with-hole mesh: 221 nodes, 192 QUAD4 cells; a = 0.1 m, W = 1.0 m, σ_∞ = 10.0 MPa
E = 210 GPa, ν = 0.3, a/W = 0.10

Static solve

res = m.solve()

Extract σ_yy at the hole edge — the Kirsch peak

sigma_yy_max = problem.extract(m, res, "kt_at_hole")
sigma_yy_kirsch = 3.0 * problem.sigma_infty
kt_computed = sigma_yy_max / problem.sigma_infty
err_pct = (sigma_yy_max - sigma_yy_kirsch) / sigma_yy_kirsch * 100.0

print()
print(f"σ_yy at hole edge   computed = {sigma_yy_max / 1e6:8.3f} MPa")
print(f"σ_yy at hole edge   Kirsch   = {sigma_yy_kirsch / 1e6:8.3f} MPa  (3 · σ_∞)")
print(f"K_t computed                 = {kt_computed:8.4f}")
print("K_t Kirsch (infinite plate)  = 3.0000")
print("K_t Peterson (a/W = 0.1)     ≈ 3.04 (Howland 1930 finite-width correction)")
print(f"relative error               = {err_pct:+.2f} %")

# 10% tolerance — Kirsch's reference is for an infinite plate; the
# finite-width Howland correction adds ~1 % at a/W = 0.1, and a
# QUAD4 mesh at 16 × 12 hits the gradient with ~3-4 % discretisation
# error on top.
assert abs(err_pct) < 10.0, f"K_t deviation {err_pct:.2f}% exceeds 10%"

σ_yy at hole edge   computed =   31.418 MPa
σ_yy at hole edge   Kirsch   =   30.000 MPa  (3 · σ_∞)
K_t computed                 =   3.1418
K_t Kirsch (infinite plate)  = 3.0000
K_t Peterson (a/W = 0.1)     ≈ 3.04 (Howland 1930 finite-width correction)
relative error               = +4.73 %

Render the σ_yy field on the deformed mesh

grid = m.grid.copy()
u = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 2)
disp3 = np.zeros((grid.n_points, 3), dtype=np.float64)
disp3[:, :2] = u
grid.point_data["displacement"] = disp3
grid.point_data["|u|"] = np.linalg.norm(u, axis=1)

# σ_yy approximation across the mesh: derive from the displacement
# field via a finite-difference of u_y in y at every interior node.
# Adequate for a contour render; the validation problem's `extract`
# method computes σ_yy(hole edge) more carefully via the plane-stress
# C-matrix at a single key point.
pts = np.asarray(grid.points)
sigma_yy_field = np.zeros(grid.n_points, dtype=np.float64)
# Use the same scaling factor from the analytical field at the hole
# edge so the colour bar's range covers [0, K_t · σ_∞] cleanly.
sigma_yy_field[:] = problem.sigma_infty  # stays at σ_∞ far from the hole
# Spike σ_yy near the hole edge (r ≈ a) to the analytical
# 3·σ_∞ value at θ = π/2.  This is for visualisation only — the
# rigorous value comes from problem.extract.
r = np.linalg.norm(pts[:, :2], axis=1)
theta = np.arctan2(pts[:, 1], pts[:, 0])
# Kirsch closed-form σ_yy(r, θ) on the y-axis (cos 2θ = -1 at θ = π/2):
near_hole = r > 0
ratio = problem.a**2 / r[near_hole] ** 2
ratio4 = problem.a**4 / r[near_hole] ** 4
sigma_yy_field[near_hole] = problem.sigma_infty * (
    1.0 + 0.5 * ratio + (1.5 * ratio4 - 0.5 * ratio) * np.cos(2.0 * theta[near_hole])
)
grid.point_data["σ_yy [Pa]"] = sigma_yy_field

warped = grid.warp_by_vector("displacement", factor=200.0)

plotter = pv.Plotter(off_screen=True, window_size=(720, 480))
plotter.add_mesh(
    grid,
    style="wireframe",
    color="grey",
    opacity=0.35,
    line_width=1,
    label="undeformed",
)
plotter.add_mesh(
    warped,
    scalars="σ_yy [Pa]",
    cmap="plasma",
    show_edges=True,
    scalar_bar_args={"title": "σ_yy [Pa]  (deformed ×200, Kirsch closed form)"},
)
# Mark the Kirsch peak point at (0, a) — the locus of K_t = 3.
plotter.add_points(
    np.array([[0.0, problem.a, 0.0]]),
    render_points_as_spheres=True,
    point_size=18,
    color="#d62728",
    label=f"K_t = {kt_computed:.3f} at (0, a)",
)
plotter.view_xy()
plotter.camera.zoom(1.05)
plotter.show()