RasmahRasmah

Plasticity of a Steel Tensile Bar (J2)

using Rasmah
using LinearAlgebra

Overview

Elasticity assumes the material comes back; plasticity is the physics of what happens when it does not. A ductile metal loaded past its yield point accumulates permanent strain, and a good analysis must predict both how much it stretches and where the plastic zone spreads. The standard model is J2 plasticity: yield starts when the von Mises (equivalent) stress reaches the yield strength σy, and hardening makes the yield surface grow as plastic strain accumulates.

This tutorial pulls a steel tensile bar — the 100 × 20 × 4 mm coupon a materials lab tests in a universal testing machine — past its yield point. Plasticity() solves the small-strain J2 problem; the load is applied in one step and the solution is found by the secant-modulus iteration, which is path-correct for monotonic loading (exactly the tensile test). The bar has an exact answer: away from the grips every cross-section follows the same uniaxial stress–strain curve, so the strain in the gauge section — the central part of the bar, where the field is uniform — must equal the strain read off the J2 hardening law. That check closes the tutorial.

The workflow is the canonical one, on a box part so face selection happens on the exact BRep before meshing:

shape → boundary sets → mesh → boundary conditions → solve

The uniaxial J2 curve

Under uniaxial tension the von Mises stress equals the axial stress, so the 1-D response is piecewise: elastic σ = E ε up to σy, then

\[ \varepsilon = \frac{\sigma_y}{E} + \frac{\sigma - \sigma_y}{H} \qquad(\sigma > \sigma_y), \]

with H the hardening modulus (the slope of the plastic branch). Below yield the bar stretches by σL/E; above it, the second term is permanent strain. The steel here has σy = 250 MPa and H = 2 GPa — a mild work-hardening structural steel.

Geometry, face selection, and meshing

The coupon is a flat box, gripped at one end and pulled at the other:

part = box(0.100, 0.020, 0.004)                       # 100 × 20 × 4 mm bar (m)

b = brep(part)
left  = faces_in_box(b, [-0.051, -0.011, -0.003], [-0.049, 0.011, 0.003])   # gripped end
right = faces_in_box(b, [ 0.049, -0.011, -0.003], [ 0.051, 0.011, 0.003])   # pulled end
mesh = tetrahedralize(b, 10)
TetMesh{Matrix{Float64}, Matrix{Int64}, Rasmah.BRepFaceLabels}([-0.05 -0.04 … 0.039999999999999994 0.05; -0.01 -0.01 … 0.01 0.01; -0.002 -0.002 … 0.002 0.002], [1 1 … 21 21; 2 13 … 54 55; 13 12 … 55 22; 46 46 … 66 66], Rasmah.BRepFaceLabels(FaceSet[FaceSet([2 5 … 3 6; 13 6 … 4 7; 14 17 … 15 18]), FaceSet([45 45 … 54 38; 46 56 … 55 49; 57 57 … 66 50]), FaceSet([3 2 … 4 6; 36 35 … 5 39; 37 36 … 38 40]), FaceSet([27 23 … 32 28; 28 24 … 33 61; 61 57 … 66 62]), FaceSet([12 12 1 1; 23 45 12 34; 56 56 45 45]), FaceSet([22 11 11 22; 33 44 22 55; 66 55 55 66])]))
vtk(b)

Elastic first

Pull with 100 MPa — below yield — and the bar should follow Hooke exactly:

u_el = solve(Plasticity(), mesh, steel;
             fixed = left => [0.0, 0.0, 0.0],
             loads = right => [100e6, 0.0, 0.0])

E = youngs_modulus(steel)
println("FEM stretch: ", round(maximum(ustrip.(u_el[1:3:end])); digits = 6), " m")
println("hand σL/E:   ", round(100e6 * 0.1 / E; digits = 6), " m")
FEM stretch: 5.0e-5 m
hand σL/E:   5.0e-5 m

The two agree — at 100 MPa the model is pure elasticity, and unloading would return the bar to its original length.

Pulling past yield

Now pull with 300 MPa. The bar enters the plastic branch: the total strain is the elastic part σy/E plus the plastic part (σ − σy)/H:

σy = yield_strength(steel)
H  = hardening_modulus(steel)
ε  = σy / E + (300e6 - σy) / H

u = solve(Plasticity(), mesh, steel;
          fixed = left => [0.0, 0.0, 0.0],
          loads = right => [300e6, 0.0, 0.0])

println("permanent strain: ", round((ε - 300e6 / E) * 100; digits = 2), " %")
permanent strain: 2.48 %

The whole coupon has yielded: about 2.5 % of the strain is permanent — cut the load and the bar stays longer.

umag = [norm(ustrip.(u)[(3i-2):(3i)]) for i in 1:size(mesh.nodes, 2)]
vtk(mesh; field = umag, fieldname = "|u| (plastic)")

The gauge-section check

Like a real extensometer, measure the strain on the centreline of the gauge section — between x = 0.01 and x = 0.04, far from both grips, where the stress state is the uniform uniaxial one. There the stretch must match the J2 curve exactly:

u = ustrip(u)
i1 = argmin([abs(mesh.nodes[1, i] - 0.01) + abs(mesh.nodes[2, i]) + abs(mesh.nodes[3, i]) for i in axes(mesh.nodes, 2)])
i2 = argmin([abs(mesh.nodes[1, i] - 0.04) + abs(mesh.nodes[2, i]) + abs(mesh.nodes[3, i]) for i in axes(mesh.nodes, 2)])

println("FEM gauge strain: ", round((u[3i2 - 2] - u[3i1 - 2]) / 0.03; digits = 6))
println("hand J2 curve:   ", round(ε; digits = 6))
FEM gauge strain: 0.026913
hand J2 curve:   0.02625

The gauge strain matches the 1-D J2 curve to a couple of percent. The small difference is the grip effect: clamping the whole end face also restrains the lateral Poisson contraction there, so the corners yield first and the field bends slightly near the grips — the same reason a real extensometer is mounted on the gauge section, not on the shoulders. (The maximum displacement, at the gripped corners, runs well ahead of the gauge value — a stress concentration, not a load-bearing number.)

Notes

  • yield_strength and hardening_modulus are read from the material (with stiffness for the elastic tensor C).
  • The solver is the secant-modulus iteration: each iteration assembles a linear problem with the current per-element secant stiffness, so the fixed point is the J2 solution for monotonic (proportional) loading — the standard tensile-test condition.
  • The model is small-strain J2 with linear isotropic hardening; for kinematic (Bauschinger) hardening or pressure-dependent yield, see the dedicated constitutive tutorials.

Next steps