# Von mises plasticity displacement control problem

I am trying to convert the plasticity tutorial given in Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation to displacement control loading on a square plate. I am able to get the solutions only when displacements are really small( Also when I introduce a crack, this gets worse). Can someone help me to write a better code?

from dolfin import *
from matplotlib import pyplot
import numpy as np
import warnings

# material parameters
E = Constant(180e3)
nu = Constant(0.28)
lmbda = Constant(89633)
mu = Constant(70300)
sig0 = Constant(443)  # yield strength
H = Constant(300)  # hardening modulus

mesh = UnitSquareMesh(200,200)

top = CompiledSubDomain("near(x[1], 1.0) && on_boundary")
bot = CompiledSubDomain("near(x[1], 0.0) && on_boundary")

deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
W = FunctionSpace(mesh, We)
W0 = FunctionSpace(mesh, W0e)

load = Expression("t", t = 0.0, degree=1)

bcbot= DirichletBC(V, Constant((0.0,0.0)), bot)
bc_u = [bcbot,bctop]

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries, 1)

ds = Measure("ds")(subdomain_data= boundaries)
n = FacetNormal(mesh)

sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)

def eps(v):
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def sigma(eps_el):
return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])

ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s
return as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp

def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)

# --------------------------------------------
# Global problem and Newton-Raphson procedure
# --------------------------------------------

a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm

def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return

u_r = 0.03
Nitermax, tol = 100, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
print("time step", Nincr)
A, Res = assemble_system(a_Newton, res, bc_u)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a_Newton, res, bc_u)
nRes = Res.norm("l2")
print("    Residual:", nRes)
for bc in bc_u:
bc.homogenize()
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)

is it just the problem of choosing the right time step and mesh, or am I making any serious mistake?

i am attaching the post-processing part below

Vs = FunctionSpace(mesh, "DG", 1)
PP = FunctionSpace(mesh, "DG", 0)

# VON MISES EQUIVALENT STRESS
sig_n1 = as_3D_tensor(sig)
s1 = dev(sig_n1)
sig_eq1 = local_project(sqrt(3/2.*inner(s1, s1)), Vs)
# plastic strain
ep_eq = local_project(sqrt(2/3)*p, PP)