Pinched ring — Castigliano diametrical deflection

Classical thin-circular-ring elasticity benchmark: a ring of mean radius \(R\) and rectangular cross-section \(b \times t\) (with \(b \ll R\) and \(t \ll R\)) is loaded by two equal opposing point loads \(P\) along a diameter. Castigliano’s theorem applied to the curved-beam strain energy gives the diametrical deflection along the loading axis (Timoshenko & Young 1968 §79; Roark Table 9 case 1):

(1)\[\delta_{\mathrm{diam}} \;=\; \frac{P\, R^{3}}{E\, I}\, \left(\frac{\pi}{4} - \frac{2}{\pi}\right),\]

with \(I = b\, t^{3} / 12\) the curved-beam bending inertia about the cross-section’s neutral axis perpendicular to the ring plane.

A quarter-symmetry FE model exploits the two orthogonal mirror planes of the loaded ring; the loaded point’s in-plane displacement equals \(\delta_{\mathrm{diam}} / 2\) since each half of the ring contributes one half of the diametrical motion.

Implementation

Drives the existing PinchedRing problem class. Quarter-arc of a thin ring discretised with 40 BEAM2 elements; the loaded end has a single \(P\)-direction force; the symmetry-plane end has the mirror-image rotation / out-of-plane DOFs pinned.

The classical Castigliano derivation uses inextensible-ring approximations (no axial stretching contribution) — at the default geometry \(R = 0.1\), \(b \times t = 0.01 \times 0.005\), the FE solution recovers the closed form to within 0.1 %, well below the 5 % tolerance the validation suite specifies.

References

• Timoshenko, S. P. and Young, D. H. (1968) Elements of Strength of Materials, 5th ed., Van Nostrand, §79 (Castigliano on a circular ring).

• Roark, R. J. and Young, W. C. (1989) Roark’s Formulas for Stress and Strain, 6th ed., McGraw-Hill, Table 9 case 1 (closed circular ring under two opposing point loads).

• Bickford, W. B. (1968) Advanced Mechanics of Materials, Addison-Wesley, §6.4 (curved-beam strain energy).

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §16.5 (curved-beam FE).

Cross-references (problem-id only):

• Abaqus AVM 1.5.x Pinched ring.

from __future__ import annotations

import math

import numpy as np
import pyvista as pv

from femorph_solver.validation.problems import PinchedRing

Build the model from the validation problem class

problem = PinchedRing()
m = problem.build_model()

# Quarter-arc geometry: 40 BEAM2 elements along π/2 of the ring.
print(
    f"pinched ring quarter-arc: {m.grid.n_points} nodes, {m.grid.n_cells} BEAM2 cells; "
    f"R = {problem.R} m, b × t = {problem.width} × {problem.height} m, P = {problem.P} N"
)
print(f"E = {problem.E / 1e9:.0f} GPa, ν = {problem.nu}")

pinched ring quarter-arc: 41 nodes, 40 BEAM2 cells; R = 0.1 m, b × t = 0.01 × 0.005 m, P = 10.0 N
E = 200 GPa, ν = 0.3

Closed-form Castigliano reference

I_section = problem.width * problem.height**3 / 12.0
delta_full = problem.P * problem.R**3 / (problem.E * I_section) * (math.pi / 4.0 - 2.0 / math.pi)
ux_quarter = delta_full / 2.0  # quarter-symmetry split

print(f"I (b·t³/12) = {I_section:.3e} m^4")
print(f"δ_diam (full ring, Timoshenko §79)   = {delta_full * 1e6:.3f} µm")
print(f"UX at loaded point (quarter model)    = {ux_quarter * 1e6:.3f} µm")

I (b·t³/12) = 1.042e-10 m^4
δ_diam (full ring, Timoshenko §79)   = 71.414 µm
UX at loaded point (quarter model)    = 35.707 µm

Static solve + extract

res = m.solve()
ux_computed = problem.extract(m, res, "loaded_point_ux")
err_pct = (ux_computed - ux_quarter) / ux_quarter * 100.0

print()
print(f"UX computed       = {ux_computed * 1e6:8.4f} µm")
print(f"UX published      = {ux_quarter * 1e6:8.4f} µm   (Castigliano)")
print(f"relative error    = {err_pct:+.4f} %")

# 0.5% tolerance — BEAM2 with 40 segments captures the curved-beam
# strain energy to within ~0.1 % on this geometry; the validation
# class permits 5 % to soak up coarser-mesh runs.
assert abs(err_pct) < 0.5, f"deflection deviation {err_pct:.4f}% exceeds 0.5 %"

UX computed       =  35.7346 µm
UX published      =  35.7068 µm   (Castigliano)
relative error    = +0.0779 %

Render the deformed quarter-arc + the analytical diametral motion

dof = m.dof_map()
disp = np.asarray(res.displacement, dtype=np.float64)
pts = np.asarray(m.grid.points)
disp_xyz = np.zeros_like(pts)
for i in range(pts.shape[0]):
    rows = np.where(dof[:, 0] == i + 1)[0]
    for r in rows:
        d_idx = int(dof[r, 1])
        if d_idx < 3:
            disp_xyz[i, d_idx] = float(disp[r])

warped = m.grid.copy()
warped.points = pts + 1.0e3 * disp_xyz  # ×1000 warp factor

plotter = pv.Plotter(off_screen=True, window_size=(640, 640))
plotter.add_mesh(m.grid, color="grey", line_width=2, opacity=0.4, label="undeformed quarter")
plotter.add_mesh(warped, color="#1f77b4", line_width=4, label="deformed (×1000 warp)")

# Mark the two characteristic loci: loaded point (P applied) at
# (R, 0) and the symmetry-plane apex (0, R).
plotter.add_points(
    np.array([[problem.R, 0.0, 0.0]]),
    render_points_as_spheres=True,
    point_size=18,
    color="#d62728",
    label=f"loaded point — UX = {ux_computed * 1e6:.3f} µm",
)
plotter.add_points(
    np.array([[0.0, problem.R, 0.0]]),
    render_points_as_spheres=True,
    point_size=14,
    color="black",
    label="symmetry plane (UY-mirror)",
)
plotter.add_legend()
plotter.view_xy()
plotter.camera.zoom(1.05)
plotter.show()