Plasticity with Isotropic Hardening#

Binder

Small-strain Plasticity

  • linear-elastic plastic material formulation with isotropic hardening

  • define a body force vector

  • extract state variables

The normal plastic strain due to a body force applied on a solid with a linear-elastic plastic material formulation with isotropic hardening with young’s modulus \(E=1\), poisson ratio \(\nu=0.3\), isotropic hardening modulus \(K=0.3\), yield stress \(\sigma_y=1\), length \(L=3\) and cross section area \(A=1\) is to be evaluated.

plasticity

First, let’s create a meshed cube out of hexahedron cells with n=(16, 6, 6) points per axis. A three-dimensional vector-valued displacement field is initiated on the numeric region.

import numpy as np
import felupe as fem

mesh = fem.Cube(b=(3, 1, 1), n=(16, 6, 6))
region = fem.RegionHexahedron(mesh)
displacement = fem.Field(region, dim=3)
field = fem.FieldContainer([displacement])

A fixed boundary condition is applied at \(x=0\).

boundaries = {"fixed": fem.dof.Boundary(displacement, fx=0)}

The material behaviour is defined through a built-in isotropic linear-elastic plastic material formulation with isotropic hardening.

Linear-elastic-plastic material formulation with linear isotropic hardening (return mapping algorithm)

  1. Given state in point x (σn, ζn=[εpn, αn]) (valid).

  2. Given strain increment dε, so that ε = εn + dε.

  3. Evaluation of the hypothetic trial state:

    dσdε = λ 1 ⊗ 1 + 2μ 1 ⊙ 1

    σ = σn + dσdε : dε

    s = dev(σ)

    εp = εpn

    α = αn

    f = ||s|| - sqrt(2/3) (σy + K α)

  4. If f ≤ 0, then elastic step:

    Set y = yn + dy, y=(σ, ζ=[εp, α]),

algorithmic consistent tangent modulus dσdε.

    Else:

    dγ = f / (2μ + 2/3 K)

n = s / ||s||

σ = σ - 2μ dγ n

εp = εpn + dγ n

α = αn + sqrt(2 / 3) dγ

Algorithmic consistent tangent modulus:

dσdε = dσdε - 2μ / (1 + K / 3μ) n ⊗ n
     - 2μ dγ / ||s|| ((2μ 1 ⊙ 1 - 1/3 1 ⊗ 1) - 2μ n ⊗ n)

umat = fem.LinearElasticPlasticIsotropicHardening(E=1, nu=0.3, sy=1, K=0.3)
solid = fem.SolidBody(umat, field)

The body force is created on the field of the solid body.

\[\delta W_{ext} = \int_v \delta \boldsymbol{u} \cdot \boldsymbol{b} ~ dv\]
bodyforce = fem.SolidBodyGravity(field)
b = fem.math.linsteps([0, 0.5], num=10, axis=0, axes=3)

Inside a Newton-Rhapson procedure, the vectors and matrices are assembled for the given items and the boundary conditions are incorporated into the equilibrium equations.

step = fem.Step(items=[solid, bodyforce], ramp={bodyforce: b}, boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate()

A view on the field-container shows the deformed mesh and the normal plastic strain in direction \(x\) due to the applied body force. The vector of all state variables, stored as a result in the solid body object, are splitted into separate variables. The plastic strain is stored as the second state variable. The mean-per-cell value of the plastic strain in direction \(x\) is exported to the view.

def screenshot(solid, field, filename):
    "Visualize the final Normal Plastic Strain in direction X."

    plastic_strain = np.split(solid.results.statevars, umat.statevars_offsets)[1]
    view = fem.View(field, cell_data={"Plastic Strain": plastic_strain.mean(-2).T})
    plotter = view.plot("Plastic Strain", component=0, off_screen=True)
    plotter.show(screenshot=filename, window_size=(800, 600), jupyter_backend="panel")

screenshot(solid, field, filename="images/bodyforce-load.png")

bodyforce-load

step = fem.Step(
    items=[solid, bodyforce], ramp={bodyforce: b[::-1]}, boundaries=boundaries
)

job = fem.Job(steps=[step])
job.evaluate()

screenshot(solid, field, filename="images/bodyforce-unload.png")

bodyforce-unload