BeamBending

A rectangular section bent past yield.

The plastic front moves in from the outer fibres, leaving an elastic core of half-depth \(c = \varepsilon_y/\kappa\). Integrating the stress over the section gives

\[\frac{M}{M_e} = \frac32 - \frac12\left(\frac{\kappa_e}{\kappa}\right)^2\]

with \(M_e = \sigma_y w h^2/6\) at first yield and \(M_p = \sigma_y w h^2/4\) once the core has vanished. Their ratio is the shape factor, exactly 3/2 for a rectangle — independent of material, size and mesh.

The ends are given a linear axial displacement, which imposes the plane-sections kinematics the closed form assumes, so this tests the constitutive response rather than beam theory. The sweep stops at \(\kappa/\kappa_e = 4\): with no hardening the curve approaches 3/2 without reaching it.

Reference

Chakrabarty, Theory of Plasticity, 3rd ed., Elsevier (2006), ch. 3.

Static Scene

Interactive Scene

• 
• 

M_e = 83333 N.mm, M_p = 125000 N.mm, shape factor = 1.500


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


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


======= elastic problem at iteration 2 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 4.952211528844e-03
At Newton iteration 3 relative norm is 2.371505112818e-05
At Newton iteration 4 relative norm is 5.752090419208e-10
At Newton iteration 5 relative norm is 1.166346001912e-14


======= elastic problem at iteration 3 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 7.846123816232e-03
At Newton iteration 3 relative norm is 3.977265766034e-05
At Newton iteration 4 relative norm is 1.917222017193e-08
At Newton iteration 5 relative norm is 2.337201284208e-14


======= elastic problem at iteration 4 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.228954348053e-02
At Newton iteration 3 relative norm is 1.593044545820e-03
At Newton iteration 4 relative norm is 1.254434127300e-05
At Newton iteration 5 relative norm is 1.628854576204e-09
At Newton iteration 6 relative norm is 3.675075308906e-14


======= elastic problem at iteration 5 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 8.241480261814e-03
At Newton iteration 3 relative norm is 2.351620835505e-05
At Newton iteration 4 relative norm is 2.017522115952e-10
At Newton iteration 5 relative norm is 4.693948764580e-14


======= elastic problem at iteration 6 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 2.215693215890e-03
At Newton iteration 3 relative norm is 1.958380994561e-06
At Newton iteration 4 relative norm is 2.146014267670e-12


======= elastic problem at iteration 7 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 3.349160557800e-04
At Newton iteration 3 relative norm is 5.845943724850e-09
At Newton iteration 4 relative norm is 8.326913701400e-14


======= elastic problem at iteration 8 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 1.117056869328e-03
At Newton iteration 3 relative norm is 1.642645443014e-04
At Newton iteration 4 relative norm is 2.133616243945e-06
At Newton iteration 5 relative norm is 3.092144534243e-10
At Newton iteration 6 relative norm is 1.037882058749e-13


======= elastic problem at iteration 9 =======
At Newton iteration 1 relative norm is 1.000000000000e+00
At Newton iteration 2 relative norm is 8.709189750462e-03
At Newton iteration 3 relative norm is 1.902986569058e-04
At Newton iteration 4 relative norm is 8.905060705193e-08
At Newton iteration 5 relative norm is 1.630443872220e-13
worst error over the sweep: 1.62 %
M/M_e reaches 1.4531, approaching the shape factor 1.5

Generate movie 01/10 (10.00 %) 1.10 s
Generate movie 02/10 (20.00 %) 868.56 ms
Generate movie 03/10 (30.00 %) 753.18 ms
Generate movie 04/10 (40.00 %) 640.22 ms
Generate movie 05/10 (50.00 %) 535.33 ms
Generate movie 06/10 (60.00 %) 435.31 ms
Generate movie 07/10 (70.00 %) 326.01 ms
Generate movie 08/10 (80.00 %) 220.83 ms
Generate movie 09/10 (90.00 %) 110.59 ms
Generate movie 10/10 (100.00 %) 0.00 µs

import numpy as np
from scipy.integrate import simpson

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

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

L, h, w = 80.0, 20.0, 5.0  # mm — length, depth, thickness
E, v = 210000.0, 0.3  # MPa
sigma_y = 250.0  # MPa

M_e = sigma_y * w * h**2 / 6  # first yield
M_p = sigma_y * w * h**2 / 4  # fully plastic
kappa_e = 2 * sigma_y / (E * h)  # curvature at first yield

print(f"M_e = {M_e:.0f} N.mm, M_p = {M_p:.0f} N.mm, shape factor = {M_p / M_e:.3f}")


def Exact(ratio):
    """M/M_e as a function of kappa/kappa_e; linear until the outer fibre yields."""
    return np.where(ratio <= 1.0, ratio, 1.5 - 0.5 / np.maximum(ratio, 1e-12) ** 2)


# ----------------------------------------------
# Mesh
# ----------------------------------------------
domain = Domain((0, -h / 2), (L, h / 2), h / 8)
mesh = domain.Mesh_2D(
    [],
    ElemType.QUAD9,
    additionalLines=[Line((0, 0), (L, 0))],
    isOrganised=True,
)
# mesh = domain.Mesh_2D([], ElemType.QUAD8, isOrganised=True)

nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
nodesXL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
nodesMid = mesh.Nodes_Conditions(lambda x, y, z: np.isclose(x, L / 2))
origin = mesh.Nodes_Conditions(lambda x, y, z: (x == 0) & (np.abs(y) < 1e-9))

yMid = mesh.coord[nodesMid, 1]
order = np.argsort(yMid)
yMid = yMid[order]

# ----------------------------------------------
# Simulation
# ----------------------------------------------
material = Models.InElastic.Behavior(
    2,
    Isotropic(3, E=E, v=v),
    yieldSurface=Models.InElastic.Yield.VonMises(sigma_y),  # perfectly plastic
    thickness=w,
    planeStress=True,
)
simu = Simulations.InElastic(mesh, material)

# Bend it, curvature by curvature
ratios = np.arange(0.4, 4.01, 0.4)
moments, profiles = [], {}

for ratio in ratios:
    kappa = ratio * kappa_e
    simu.Bc_Init()
    simu.add_dirichlet(nodesX0, [0], ["x"])
    simu.add_dirichlet(origin, [0], ["y"])
    # plane sections stay plane: u_x = -kappa * L * y at the far end
    simu.add_dirichlet(nodesXL, [lambda x, y, z: -kappa * L * y], ["x"])
    simu.Solve()
    simu.Save_Iter()

    sig_xx = simu.Result("Sxx")[nodesMid][order]
    # Simpson, since sig_xx * y is quadratic in y while the section is elastic
    moments.append(abs(simpson(sig_xx * yMid, x=yMid) * w))
    profiles[ratio] = sig_xx

moments = np.array(moments)
exact = Exact(ratios)

errors = 100 * np.abs(moments / M_e - exact) / exact

shapeFactor = moments[-1] / M_e
print(f"worst error over the sweep: {errors.max():.2f} %")
print(f"M/M_e reaches {shapeFactor:.4f}, approaching the shape factor 1.5")

# only the shape factor is asserted: it is a pure number, independent of material, size and
# mesh, where the error against the closed form moves with any mesh change
assert 1.4 < shapeFactor < 1.5, "the shape factor is not approached"

# ----------------------------------------------
# Results
# ----------------------------------------------
PyVista.Plot_BoundaryConditions(simu).show()

# the two curves lie on top of each other, so the error gets a panel of its own
fig, (ax, axErr) = Matplotlib.plt.subplots(
    2, 1, sharex=True, figsize=(6.4, 6.0), gridspec_kw={"height_ratios": [3, 1]}
)

fine = np.linspace(0.05, ratios[-1], 400)
ax.axvspan(0, 1, color="0.93", zorder=0)
ax.text(0.5, 0.07, "elastic", ha="center", fontsize=8, color="0.4")
ax.plot(fine, Exact(fine), "k-", lw=1.2, label="Chakrabarty")
ax.plot(ratios, moments / M_e, "o", ms=6, mfc="none", mew=1.3, label="EasyFEA")
ax.axhline(1.5, ls="--", c="r", lw=0.9)
ax.text(
    ratios[-1], 1.505, "shape factor $M_p/M_e = 3/2$  ", c="r", fontsize=8, ha="right"
)
ax.set_ylabel("$M/M_e$")
ax.set_ylim(0, 1.62)
ax.set_xlim(0, ratios[-1] + 0.1)
ax.set_title("Moment-curvature of a rectangular section")
ax.legend(fontsize=8, loc="lower right")
ax.grid(alpha=0.3)

axErr.axvspan(0, 1, color="0.93", zorder=0)
axErr.plot(ratios, errors, "o-", ms=4, lw=1)
axErr.set_xlabel(r"$\kappa/\kappa_e$")
axErr.set_ylabel("error [%]")
axErr.set_ylim(0, max(1.05 * errors.max(), 0.5))
axErr.grid(alpha=0.3)
fig.align_ylabels()

# the plastic front eating into the section
ax = Matplotlib.Init_Axes()
yFine = np.linspace(-h / 2, h / 2, 400)
for i, ratio in enumerate((0.8, 1.6, 2.4, 4.0)):
    idx = int(np.argmin(np.abs(ratios - ratio)))
    ax.plot(
        profiles[ratios[idx]], yMid, label=rf"$\kappa/\kappa_e$ = {ratios[idx]:.1f}"
    )
    # exact: linear inside the elastic core of half-depth c, saturated at sigma_y outside.
    # the imposed u_x = -kappa L y puts the top fibre in compression, hence the sign
    core = h / (2 * ratios[idx])
    ax.plot(
        -sigma_y * np.clip(yFine / core, -1, 1),
        yFine,
        "k--",
        lw=0.8,
        label="exact" if i == 0 else None,
    )
ax.axvline(sigma_y, ls=":", c="k", lw=0.8)
ax.axvline(-sigma_y, ls=":", c="k", lw=0.8)
ax.set_xlabel(r"$\sigma_{xx}$ [MPa]")
ax.set_ylabel("$y$ [mm]")
ax.set_title("Stress through the section: the elastic core shrinks")
ax.legend()
ax.grid(alpha=0.3)

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

Matplotlib.plt.show()