ThickCylinder

The classic elastic-plastic benchmark: a thick-walled cylinder under internal pressure, with three independent analytical landmarks on one problem.

For an elastic-perfectly-plastic material in plane strain the von Mises condition reduces to \(\sigma_\theta - \sigma_r = Y\) with \(Y = 2\sigma_y/\sqrt3\), and equilibrium gives the pressure putting the elastic-plastic boundary at radius c:

\[p(c) = \frac{Y}{2}\left[\,2\ln\frac{c}{a} + 1 - \frac{c^2}{b^2}\right]\]

Below \(p_e = \frac{Y}{2}(1 - a^2/b^2)\) the wall is elastic and Lamé holds. Between the two the stresses follow from equilibrium and the yield condition alone, so no elastic constant enters the plastic zone. At \(p_{lim} = Y\ln(b/a)\) the wall is fully plastic and the cylinder collapses.

The elastic-plastic boundary is a kink that no mesh resolves exactly, so the error converges only first order even with quadratic elements.

References

Hill, The Mathematical Theory of Plasticity, Oxford (1950), ch. V.

Bleyer, Elasto-plastic analysis of a 2D von Mises material, Computational Mechanics Numerical Tours with FEniCSx — the same cylinder, hardening where this one is perfectly plastic.

Static Scene

Interactive Scene

Static Scene

Interactive Scene

• 
• 

yield starts at p_e   =  108.25 MPa
fully plastic at p_lim=  200.09 MPa
applied      p      =  160.00 MPa -> plastic front at c = 130.71 mm


======= elastic problem at iteration 0 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.875144294267e-13


======= elastic problem at iteration 1 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 5.007716556554e-13


======= elastic problem at iteration 2 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 7.247492722967e-13


======= elastic problem at iteration 3 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 9.458592762635e-13


======= elastic problem at iteration 4 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.160138161556e-12


======= elastic problem at iteration 5 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.407567456597e-12


======= elastic problem at iteration 6 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.656234243374e-12


======= elastic problem at iteration 7 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 6.709404831177e-01
At Newton iteration 3 relative norm is 1.762061810711e-03
At Newton iteration 4 relative norm is 2.579457230361e-08
At Newton iteration 5 relative norm is 2.001219334515e-12


======= elastic problem at iteration 8 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 8.871841952830e-01
At Newton iteration 3 relative norm is 1.473671457143e-02
At Newton iteration 4 relative norm is 1.516873750904e-05
At Newton iteration 5 relative norm is 1.725181996731e-11


======= elastic problem at iteration 9 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.511415350518e+00
At Newton iteration 3 relative norm is 8.206972642276e-02
At Newton iteration 4 relative norm is 1.699234160099e-03
At Newton iteration 5 relative norm is 1.802452071315e-07
At Newton iteration 6 relative norm is 2.486064777166e-12


======= elastic problem at iteration 10 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.697073911286e+00
At Newton iteration 3 relative norm is 8.396749071630e-02
At Newton iteration 4 relative norm is 3.712583517307e-03
At Newton iteration 5 relative norm is 3.075750257560e-07
At Newton iteration 6 relative norm is 2.692623498406e-12


======= elastic problem at iteration 11 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.000226056590e+00
At Newton iteration 3 relative norm is 1.861389512284e-01
At Newton iteration 4 relative norm is 6.829637944726e-03
At Newton iteration 5 relative norm is 3.209380850956e-06
At Newton iteration 6 relative norm is 5.053819419793e-12


======= elastic problem at iteration 12 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.216877265179e+00
At Newton iteration 3 relative norm is 2.278110492143e-01
At Newton iteration 4 relative norm is 1.642969211136e-02
At Newton iteration 5 relative norm is 2.946412172108e-06
At Newton iteration 6 relative norm is 4.190338506385e-12


======= elastic problem at iteration 13 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.324725740337e+00
At Newton iteration 3 relative norm is 3.841760017385e-01
At Newton iteration 4 relative norm is 6.375027465279e-02
At Newton iteration 5 relative norm is 6.769390813836e-05
At Newton iteration 6 relative norm is 1.732502954327e-10


======= elastic problem at iteration 14 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.970320111005e+00
At Newton iteration 3 relative norm is 4.290808756303e-01
At Newton iteration 4 relative norm is 4.861347111367e-02
At Newton iteration 5 relative norm is 8.509775836130e-06
At Newton iteration 6 relative norm is 6.881614950666e-12


======= elastic problem at iteration 15 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.680547377544e+00
At Newton iteration 3 relative norm is 5.044231477387e-01
At Newton iteration 4 relative norm is 5.898269567130e-02
At Newton iteration 5 relative norm is 1.820417515046e-05
At Newton iteration 6 relative norm is 1.694891760301e-11


======= elastic problem at iteration 16 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.476155185380e+00
At Newton iteration 3 relative norm is 8.602336230811e-01
At Newton iteration 4 relative norm is 2.168754676369e-01
At Newton iteration 5 relative norm is 4.528916616383e-04
At Newton iteration 6 relative norm is 6.038376390773e-10


======= elastic problem at iteration 17 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.516024959295e+00
At Newton iteration 3 relative norm is 6.224310049497e-01
At Newton iteration 4 relative norm is 2.621836194463e-01
At Newton iteration 5 relative norm is 1.445262236300e-02
At Newton iteration 6 relative norm is 3.887661338086e-06
At Newton iteration 7 relative norm is 5.605731448186e-12


======= elastic problem at iteration 18 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 3.054898529550e+00
At Newton iteration 3 relative norm is 7.246032488391e-01
At Newton iteration 4 relative norm is 2.190613164383e-01
At Newton iteration 5 relative norm is 1.584180201424e-02
At Newton iteration 6 relative norm is 1.233201768484e-05
At Newton iteration 7 relative norm is 2.708432564889e-11


======= elastic problem at iteration 19 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.960792874602e+00
At Newton iteration 3 relative norm is 8.535522490082e-01
At Newton iteration 4 relative norm is 2.186122473292e-01
At Newton iteration 5 relative norm is 2.350205212945e-02
At Newton iteration 6 relative norm is 2.683545012217e-05
At Newton iteration 7 relative norm is 1.271237931589e-10


======= elastic problem at iteration 20 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.840297077445e+00
At Newton iteration 3 relative norm is 9.371569598673e-01
At Newton iteration 4 relative norm is 2.431086264328e-01
At Newton iteration 5 relative norm is 2.267975558954e-02
At Newton iteration 6 relative norm is 3.090766185058e-05
At Newton iteration 7 relative norm is 1.826097607575e-10


======= elastic problem at iteration 21 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 3.267334984616e+00
At Newton iteration 3 relative norm is 1.102864873286e+00
At Newton iteration 4 relative norm is 2.506806597952e-01
At Newton iteration 5 relative norm is 2.018497239108e-02
At Newton iteration 6 relative norm is 5.575288865454e-05
At Newton iteration 7 relative norm is 5.230469224443e-10


======= elastic problem at iteration 22 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 3.343425191189e+00
At Newton iteration 3 relative norm is 1.369502140561e+00
At Newton iteration 4 relative norm is 4.448018217609e-01
At Newton iteration 5 relative norm is 5.084370777743e-02
At Newton iteration 6 relative norm is 1.028633488190e-03
At Newton iteration 7 relative norm is 4.687670848969e-09


======= elastic problem at iteration 23 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 4.142712196158e+00
At Newton iteration 3 relative norm is 1.788098985347e+00
At Newton iteration 4 relative norm is 3.612909605364e-01
At Newton iteration 5 relative norm is 6.244702822672e-02
At Newton iteration 6 relative norm is 8.359080661240e-03
At Newton iteration 7 relative norm is 3.072047124869e-06
At Newton iteration 8 relative norm is 8.418617166959e-12


======= elastic problem at iteration 24 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 6.112836476475e+00
At Newton iteration 3 relative norm is 2.618066063215e+00
At Newton iteration 4 relative norm is 6.827716340465e-01
At Newton iteration 5 relative norm is 7.800320757878e-02
At Newton iteration 6 relative norm is 3.725871416291e-04
At Newton iteration 7 relative norm is 7.240085130174e-09


======= elastic problem at iteration 25 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 5.187380937323e+00
At Newton iteration 3 relative norm is 4.162471295846e+00
At Newton iteration 4 relative norm is 2.999823746071e+00
At Newton iteration 5 relative norm is 1.579742963819e+00
At Newton iteration 6 relative norm is 5.527243437268e-01
At Newton iteration 7 relative norm is 4.767750821698e-01
At Newton iteration 8 relative norm is 4.531106318509e-01
At Newton iteration 9 relative norm is 3.480431703165e-01
At Newton iteration 10 relative norm is 2.449549125872e-01
At Newton iteration 11 relative norm is 4.612958193308e-02
At Newton iteration 12 relative norm is 9.601061785321e-04
At Newton iteration 13 relative norm is 3.332349794160e-07
At Newton iteration 14 relative norm is 1.368930113419e-11

at p = 160 MPa on 1493 elements: max |error| vs Hill = 7.42 MPa = 2.57 % of Y
plastic front at r = 131.2 mm, Hill says 130.7, element size 6.2

Generate movie 01/26 (3.85 %) 3.31 s
Generate movie 02/26 (7.69 %) 2.94 s
Generate movie 03/26 (11.54 %) 2.87 s
Generate movie 04/26 (15.38 %) 2.72 s
Generate movie 05/26 (19.23 %) 2.57 s
Generate movie 06/26 (23.08 %) 2.49 s
Generate movie 07/26 (26.92 %) 2.38 s
Generate movie 08/26 (30.77 %) 2.28 s
Generate movie 09/26 (34.62 %) 2.13 s
Generate movie 10/26 (38.46 %) 2.02 s
Generate movie 11/26 (42.31 %) 2.01 s
Generate movie 12/26 (46.15 %) 1.80 s
Generate movie 13/26 (50.00 %) 1.63 s
Generate movie 14/26 (53.85 %) 1.51 s
Generate movie 15/26 (57.69 %) 1.38 s
Generate movie 16/26 (61.54 %) 1.25 s
Generate movie 17/26 (65.38 %) 1.13 s
Generate movie 18/26 (69.23 %) 1.01 s
Generate movie 19/26 (73.08 %) 883.50 ms
Generate movie 20/26 (76.92 %) 753.81 ms
Generate movie 21/26 (80.77 %) 622.34 ms
Generate movie 22/26 (84.62 %) 501.04 ms
Generate movie 23/26 (88.46 %) 378.18 ms
Generate movie 24/26 (92.31 %) 257.06 ms
Generate movie 25/26 (96.15 %) 127.58 ms
Generate movie 26/26 (100.00 %) 0.00 µs

import numpy as np
from scipy.optimize import brentq

from EasyFEA import Folder, ElemType, Models, Simulations, PyVista, Matplotlib
from EasyFEA.Geoms import CircleArc, Contour, Line, Circle
from EasyFEA.Models.Elastic._laws import Isotropic

# ----------------------------------------------
# Configuration
# ----------------------------------------------
folder = Folder.Results_Dir()

a, b = 100.0, 200.0  # mm, inner and outer radius
E, v = 210000.0, 0.3  # MPa
sigma_y = 250.0  # MPa

Y = 2 * sigma_y / np.sqrt(3)  # plane-strain yield in sigma_theta - sigma_r
p_e = Y / 2 * (1 - a**2 / b**2)  # bore starts to yield
p_lim = Y * np.log(b / a)  # whole wall plastic: collapse
pressure = 160.0  # MPa, between the two so the front sits inside the wall

c = brentq(lambda c: Y / 2 * (2 * np.log(c / a) + 1 - c**2 / b**2) - pressure, a, b)

print(f"yield starts at p_e   = {p_e:7.2f} MPa")
print(f"fully plastic at p_lim= {p_lim:7.2f} MPa")
print(f"applied      p      = {pressure:7.2f} MPa -> plastic front at c = {c:.2f} mm")


def Exact(r):
    """Radial and hoop stress, Hill 1950 ch. V."""
    plastic = r <= c
    sig_r = np.where(
        plastic,
        -pressure + Y * np.log(r / a),
        Y * c**2 / (2 * b**2) * (1 - b**2 / r**2),
    )
    sig_t = np.where(plastic, sig_r + Y, Y * c**2 / (2 * b**2) * (1 + b**2 / r**2))
    return sig_r, sig_t


def Elastic_u(p):
    """Bore displacement, Lamé in plane strain."""
    return (1 + v) * a * p / (E * (b**2 - a**2)) * ((1 - 2 * v) * a**2 + b**2)


# ----------------------------------------------
# Mesh
# ----------------------------------------------
origin = (0, 0)
p1 = (a, 0)
p2 = (b, 0)
p3 = (0, b)
p4 = (0, a)
meshSize = (b - a) / 16  # 16 elements through the wall
contour = Contour(
    [
        Line(p1, p2, meshSize),
        CircleArc(p2, p3, center=origin, meshSize=meshSize),
        Line(p3, p4, meshSize),
        CircleArc(p4, p1, center=origin, meshSize=meshSize),
    ]
)
mesh = contour.Mesh_2D([], ElemType.TRI6)

nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
nodesY0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
bore = mesh.Nodes_Circle(Circle((0, 0), diam=2 * a))

# perfectly plastic, as Hill assumes
material = Models.InElastic.Behavior(
    2,
    Isotropic(3, E=E, v=v),
    yieldSurface=Models.InElastic.Yield.VonMises(sigma_y),
    planeStress=False,
)


# 1 & 3. One ramp, from first yield to collapse
# ----------------------------------------------
# sqrt spacing: the increments shorten as p_lim is approached, where they must. `pressure` is
# inserted rather than jumped to, since plasticity is path dependent and the wall is read there.
steps = np.sort(np.append(p_lim * 0.998 * np.linspace(0, 1, 26)[1:] ** 0.5, pressure))
iPressure = int(np.flatnonzero(steps == pressure)[0])

simu = Simulations.InElastic(mesh, material)
node_a = mesh.Nodes_Conditions(lambda x, y, z: (y == 0) & (x <= a * 1.001))

u_bore = []
for p in steps:
    simu.Bc_Init()
    simu.add_dirichlet(nodesY0, [0], ["y"])
    simu.add_dirichlet(nodesX0, [0], ["x"])
    simu.add_pressureLoad(bore, p)
    simu.Solve()
    simu.Save_Iter()
    u_bore.append(float(np.mean(simu.Result("ux")[node_a])))
u_bore = np.array(u_bore)

# ----------------------------------------------
# 2. Stresses through a partly plastic wall
# ----------------------------------------------
simu.Set_Iter(iPressure)

# sample along y = 0, where the radial direction is x, so sig_r = Sxx and sig_t = Syy
order = np.argsort(mesh.coord[nodesY0, 0])
r = mesh.coord[nodesY0, 0][order]
sig_r = simu.Result("Sxx")[nodesY0][order]
sig_t = simu.Result("Syy")[nodesY0][order]

exact_r, exact_t = Exact(r)
err = max(np.max(np.abs(sig_r - exact_r)), np.max(np.abs(sig_t - exact_t)))
print(
    f"\nat p = {pressure:.0f} MPa on {mesh.Ne} elements: "
    f"max |error| vs Hill = {err:.2f} MPa = {100 * err / Y:.2f} % of Y"
)
# the error above is a discretisation error, so no fixed bound on it says anything about the
# physics. What does: Hill puts the front at c, and a mesh of element size h cannot place it
# any closer than that.
p_r = simu.Result("p")[nodesY0][order]
front = r[p_r > 0].max()
print(
    f"plastic front at r = {front:.1f} mm, Hill says {c:.1f}, element size {meshSize:.1f}"
)
assert abs(front - c) < meshSize, "the plastic front is not where Hill puts it"

# ----------------------------------------------
# Results
# ----------------------------------------------
rr = np.linspace(a, b, 400)
exact_r, exact_t = Exact(rr)

PyVista.Plot_BoundaryConditions(simu).show()

ax = Matplotlib.Init_Axes()
ax.plot(rr / a, exact_t, "k-", lw=1, label=r"$\sigma_\theta$ exact")
ax.plot(rr / a, exact_r, "k--", lw=1, label=r"$\sigma_r$ exact")
ax.plot(r / a, sig_t, "o", ms=3, label=r"$\sigma_\theta$ FE")
ax.plot(r / a, sig_r, "s", ms=3, label=r"$\sigma_r$ FE")
ax.axvline(c / a, ls=":", c="k", lw=0.8)
ax.text(c / a, 0, "  plastic front", rotation=90, va="bottom")
ax.set_xlabel("$r/a$")
ax.set_ylabel("stress [MPa]")
ax.set_title(f"Thick cylinder at $p$ = {pressure:.0f} MPa")
ax.legend()
ax.grid(alpha=0.3)

ax = Matplotlib.Init_Axes()
ax.plot(u_bore, steps / p_lim, "o-", ms=3, lw=1, label="FE")
# the same pressures against the displacement Lamé predicts: the FE curve leaves it at p_e
ax.plot(Elastic_u(steps), steps / p_lim, "k--", lw=1, label="Lamé (elastic)")
ax.axhline(1.0, c="r", ls="--", lw=1)
ax.text(u_bore[0], 1.0, "$p_{lim} = Y \, \\ln(b/a)$", c="r", va="bottom")
ax.axhline(p_e / p_lim, c="k", ls=":", lw=0.8)
ax.text(u_bore[0], p_e / p_lim, "$p_e$", va="bottom")
ax.set_xlabel("radial displacement at the bore [mm]")
ax.set_ylabel("$p / p_{lim}$")
ax.set_title("Load-displacement to collapse")
ax.set_ylim(0, 1.1)
ax.legend()
ax.grid(alpha=0.3)

PyVista.Plot(simu, "Svm", plotMesh=True, nColors=11).show()

PyVista.Movie_simu(simu, "p", folder, "p.gif", deformFactor=10)

Matplotlib.plt.show()