Plasticity

We subject a 3D unit cube meshed with linear hexahedrons to a uniaxial stress state in the x-direction. The displacement is increased in increments and the material yields when the stress reaches a critical value. The material is modeled with a v. Mises yield criterion and a linear isotropic hardening law.

import matplotlib.pyplot as plt
import torch

from torchfem import Solid
from torchfem.materials import IsotropicPlasticity3D
from torchfem.mesh import cube_hexa

torch.set_default_dtype(torch.float64)

Material model

The cube is modeled with IsotropicPlasticity3D, a von Mises yield criterion with an associative flow rule integrated by a radial return mapping. The material documentation states the yield surface, the return mapping and the algorithmic tangent.

Hardening is linear here,

$$ \sigma_f(q) = \sigma_y + k q $$

with the yield stress $\sigma_y$, the equivalent plastic strain $q$ and the hardening modulus $k$.

E = 1000.0
nu = 0.3
sigma_y = 50.0
k = 100.0


# Hardening function
def sigma_f(q):
    return sigma_y + k * q


# Derivative of the hardening function
def sigma_f_prime(q):
    return k


# Elastic properties
material = IsotropicPlasticity3D(E, nu, sigma_f, sigma_f_prime)

# Generate cube
nodes, elements = cube_hexa(5, 5, 5)

box = Solid(nodes, elements, material)

# Set constraints
DL = 0.1
box.displacements[nodes[:, 0] == 1.0, 0] = DL
box.constraints[nodes[:, 0] == 0.0, 0] = True
box.constraints[nodes[:, 0] == 1.0, 0] = True
box.constraints[nodes[:, 1] == 0.5, 1] = True
box.constraints[nodes[:, 2] == 0.5, 2] = True

# Incremental loading
increments = torch.cat((torch.linspace(0.0, 1.0, 10), torch.linspace(1.0, 0.0, 10)))
u, f, σ, F, α = box.solve(increments=increments, return_intermediate=True)

Postprocessing and evaluation

The reference solution for the plastic region is given by solving the hardening rule $$ \sigma = \sigma_y + K q $$ for $q$ and substituting it into the elastic equation $$ \sigma = E (\varepsilon - q). $$ Solving for $\sigma$ gives $$ \sigma = \frac{kE}{k+E} (\varepsilon + \frac{\sigma_y}{k}). $$

ref_strain = [0.0, sigma_y / E, DL, DL - k / (k + E) * (DL + sigma_y / k)]
ref_stress = [0.0, sigma_y, k * E / (k + E) * (DL + sigma_y / k), 0]

# Average x components over elements
ε = 0.5 * (F.transpose(-1, -2) + F) - torch.eye(3)
strain = ε[:, :, 0, 0].mean(dim=1)
stress = σ[:, :, 0, 0].mean(dim=1)

plt.plot(ref_strain, ref_stress, "-", color="lightgray", linewidth=5, label="Reference")
plt.plot(strain, stress, ".-", label="FEM")
plt.xlabel("Strain ε_xx")
plt.ylabel("Stress σ_xx")
plt.title("Stress-strain curve for isotropic linear hardening")
plt.ylim(bottom=0)
plt.xlim(left=0)
plt.grid()
plt.legend()
plt.show()