PDDO-MultiPhysics Examples Solid mechanics

J2 plasticity: loading and unloading a bar

Small-strain J2 plasticity with linear isotropic hardening: a bar loaded to three times its yield strain and unloaded, with the return mapping written in the problem file.

plasticity_j2_bar.yaml 275 nodes backward-euler, 23 s

1The Problem

A bar in uniaxial tension is stretched to three times the yield strain and then brought back to zero length, which drives it through elastic loading, plastic flow with hardening, elastic unloading and reverse yielding. In uniaxial stress the exact response is known, so the computed stress can be compared with it at every step.

$$\boldsymbol{\sigma}^{\text{tr}} = \mathbb{C} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p_n), \qquad q^{\text{tr}} = \sqrt{\tfrac32\,\mathbf{s}^{\text{tr}}:\mathbf{s}^{\text{tr}}}$$
$$\Delta\gamma = \frac{\langle q^{\text{tr}} - (\sigma_y + H\alpha_n)\rangle}{3G + H}, \qquad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\text{tr}} - 2G\,\Delta\gamma\,\tfrac32\,\frac{\mathbf{s}^{\text{tr}}}{q^{\text{tr}}}$$

2The Problem File

For linear hardening the radial return has a closed form, so the stress can be written as a function of the current strain and of the plastic state at the start of the step. The plastic strain and the hardening variable are fields with an update rule, applied once per converged step: eps_p: {type: symtensor, update: "eps_p + dgam*nflow"}. Because Newton differentiates the returned stress exactly, it uses the consistent tangent and converges quadratically. The intermediate q_tr and dgam are marked cache: true, which keeps the symbolic derivatives small. The load path is a table of the end displacement against time.

The complete problem file, examples/plasticity_j2_bar.yaml 
# J2 plasticity (small strain, linear isotropic hardening) - bar in uniaxial tension,
# loaded to 3x the yield strain and unloaded.
# Monolithic return mapping: for linear hardening the radial return has a closed form, so
# the stress is written as a function of the current strain and the plastic state at the
# start of the step (eps_p, alpha):
#   trial stress  sig_tr = C : (eps - eps_p),  q_tr = von Mises(sig_tr)
#   dgam = <q_tr - (sy + Hk alpha)> / (3 G + Hk),   sigma = sig_tr - 2 G dgam nflow
# Newton differentiates sigma exactly: the consistent (algorithmic) tangent, quadratic
# convergence. The internal variables eps_p, alpha are update-rule fields, updated once per
# converged step with the same closed form. q_tr and dgam are cached (evaluated pointwise):
# much smaller symbolic derivatives. (A staggered alternative: the trial stress in
# div(sigma) with solver.staggered iterations - linear convergence, see README.)
# Uniaxial analytic curve: sigma = E eps up to sigma_y, then sigma_y + E H/(E + H) (eps - eps_y).
"k">name: plasticity_j2_bar
"k">domain:
  "k">size: [1.0, 0.2, 0.2]
  "k">n: [11, 5, 5]
"k">parameters:
  "k">E: 200.0
  "k">nu: 0.3
  "k">sy: 0.2                         # yield stress (eps_y = 1e-3)
  "k">Hk: 20.0                        # isotropic hardening modulus
  "k">lambda: "E*nu/((1 + nu)*(1 - 2*nu))"
  "k">G: "E/(2*(1 + nu))"
  "k">ubar: {table: [[0, 0], [1, 3.0e-3], [2, 0]], of: t}     # end displacement (load path)
"k">fields:
  "k">u: vector
  "k">eps_p: {type: symtensor, update: "eps_p + dgam*nflow"}
  "k">alpha: {type: scalar, update: "alpha + dgam"}
"k">relations:
  "k">eps: "sym(grad(u))"
  "k">sig_tr: "lambda*tr(eps - eps_p)*I + 2*G*(eps - eps_p)"   # trial (elastic predictor)
  "k">s_tr: "dev(sig_tr)"
  "k">q_tr: {expr: "sqrt(3/2*ddot(s_tr, s_tr) + 1e-24)", cache: true}   # von Mises (smooth at 0)
  "k">dgam: {expr: "max(q_tr - (sy + Hk*alpha), 0)/(3*G + Hk)", cache: true}   # plastic multiplier
  "k">nflow: "3/2*s_tr/q_tr"                                   # flow direction
  "k">sigma: "sig_tr - 2*G*dgam*nflow"                         # return-mapped stress
"k">equations:
  "k">u: "div(sigma) = 0"
"k">boundary_conditions:
  - "k">faces: [x-]
    "k">u: {equation: "dot(sigma, n) = [0, 0, 0]"}
    "k">u_x: {dirichlet: 0}
  - "k">faces: [y-]
    "k">u: {equation: "dot(sigma, n) = [0, 0, 0]"}
    "k">u_y: {dirichlet: 0}
  - "k">faces: [z-]
    "k">u: {equation: "dot(sigma, n) = [0, 0, 0]"}
    "k">u_z: {dirichlet: 0}
  - "k">faces: [x+]
    "k">u: {equation: "dot(sigma, n) = [0, 0, 0]"}
    "k">u_x: {dirichlet: "ubar"}
  - "k">faces: all
    "k">u: {equation: "dot(sigma, n) = [0, 0, 0]"}
"k">time:
  "k">end: 2.0
  "k">dt: 0.05
  "k">scheme: backward-euler
"k">output:
  "k">every: 10
"k">monitors:
  "k">strain: "mean(eps[0, 0])"
  "k">stress: "mean(sigma[0, 0])"
  "k">eps_p_xx: "mean(eps_p_xx)"
  "k">alpha: "max(alpha)"
The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.
Figure 1. The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.
The return mapping as relations: trial stress, von Mises stress, plastic multiplier, flow direction and returned stress.
Figure 2. The return mapping as relations: trial stress, von Mises stress, plastic multiplier, flow direction and returned stress.
Internal variables as fields with update rules: the plastic strain tensor and the hardening variable.
Figure 3. Internal variables as fields with update rules: the plastic strain tensor and the hardening variable.
Material parameters and the load path, a table of the end displacement against time.
Figure 4. Material parameters and the load path, a table of the end displacement against time.

3Results

Nodes
11×5×5 = 275
PDDO order
2, horizon 2.015 Δx
Time scheme
backward-euler
Steps
40 × Δt = 0.05
Newton iterations
94
Factorizations
95
Run time
23 s

This run stored every output step for the plots: --set output.every=1 --set "output.fields=[u, eps_p, alpha]" on the command line; the problem file itself writes less often.

The stress peaks at 0.2364 at the largest stretch and reaches -0.2479 after reverse yielding at the end of unloading. Replaying the uniaxial return mapping over the same strain history, the largest difference from the computed stress is 1.2×10-9 of the peak. Newton needed 2.4 iterations per load step on average.

Mean axial stress against time: elastic, hardening, elastic unloading and reverse yielding.
Figure 5. Mean axial stress against time: elastic, hardening, elastic unloading and reverse yielding.
Mean axial strain and plastic strain against time.
Figure 6. Mean axial strain and plastic strain against time.
The hardening variable &alpha; at the largest stretch (t = 1), on a fixed colour range from 0 to 0.0025: uniform, 0.00182, as it should be in uniaxial stress.
Figure 7. The hardening variable α at the largest stretch (t = 1), on a fixed colour range from 0 to 0.0025: uniform, 0.00182, as it should be in uniaxial stress.

4Run It Yourself

The problem file ships with PDDO-MultiPhysics. Open it in the graphical interface and press Compute, or run it from the command line:

python -m pddo_mp gui examples/plasticity_j2_bar.yaml
python -m pddo_mp run examples/plasticity_j2_bar.yaml --set output.every=1 --set "output.fields=[u, eps_p, alpha]"

Results go to results/plasticity_j2_bar/ next to the problem file: a ParaView time series (.pvd), the monitors as CSV, and the final state. Any entry can be changed on the command line, for example --set time.end=1 for a run to half the end time. The user manual describes every part of the problem file and of the interface.

Have a coupled problem of your own?

Send us the physics and the geometry — we will write the problem file with you, or build the model you need.

Talk to an engineer About PDDO-MultiPhysics