Von Mises plasticity example

dolfinx
1

Good morning, I’m trying to modify the code written by @bleyerj for the Von Mises elasto-plasticity example: Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation
considering a holed rectangular plate subjected to an incremental traction load, from 0 up to
1.15 \sigma_0. The load is applied on the top boundary.
Here I show my code:


from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)

# elastic parameters
E = Constant(70e3)
nu = Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = Constant(250.)  # yield strength
Et = E/100.  # tangent modulus
H = E*Et/(E-Et)  # hardening modulus
from mshr import *
domain = Rectangle(Point(0.,0.), Point(1., 2.)) - Circle(Point(0., 0.), 0.2)
mesh = generate_mesh(domain, 50)
plot(mesh)
def bottom(x):
    return np.isclose(x[1], 0)

def top(x):
    return np.isclose(x[1], 2)

def left(x):
  return np.isclose(x[0],0)

D = mesh.topology().dim()
print(D)
neumann_domain = MeshFunction("size_t", mesh, D-1)
neumann_domain.set_all(0)
CompiledSubDomain("near(x[1], side) && on_boundary", side=2.0, tol=10e-15).mark(neumann_domain, 1)
ds = Measure("ds", subdomain_data=neumann_domain)

deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
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)

bc = [DirichletBC(V.sub(1), 0, bottom), DirichletBC(V.sub(0), 0, left )]

n = FacetNormal(mesh)
q_lim = float(sig0*1.15)
loading = Expression("q*t", q=q_lim, t=0, degree=2)

def F_ext(v):
    return loading*dot(n, v)*ds(1)
def eps(v):
    e = sym(grad(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)
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)
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
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")

Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
    loading.t = t
    A, Res = assemble_system(a_Newton, res, bc)
    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)
        nRes = Res.norm("l2")
        print("    Residual:", nRes)
        niter += 1
    u.assign(u+Du)
    p.assign(p+local_project(dp_, W0))
    sig_old.assign(sig)
    file_results.write(u, t)
    p_avg.assign(project(p, P0))
    file_results.write(p_avg, t)

I don’t know if the code is correct or not.
I get this result for the displacement correspondent to the maximum load:

IMG_8079
IMG_8079970×1277 100 KB

Could you help me understand if the results are correct?

2

@bleyerj I would also like to ask you how to plot the stress components. I have tried with

plot(sig(u)[0,0],mode=‘color’)

but it doesn’t work.

3

You should project on a tensor function space which can be plotted such as DG0 for instance.
If your code is running don’t expect someone to help you in understanding if your results are correct or not. You should post here a real implementation issue or difficulty.

4

Thank you for the advice. The code seems to work correctly, I was not very sure about the results. I have tried also another modification of your code introducing an imposed displacement, but doesn’t seem to work. I may open another discussion for it.
Regarding the stress:

sig_n1 = as_3D_tensor(sig)
s1 = dev(sig_n1)
Vs = FunctionSpace(mesh,"DG",1)
sig_eq1 = local_project(sqrt(3/2.*inner(s1,s1)), Vs)

import matplotlib.pyplot as plt
p = plot(sig_eq1, mode='color')
plt.colorbar(p)
plt.title(r"Von Mises",fontsize=26)
plt.show()

I have used this for the Von Mises stress, that I have found in another discussion.

And this seems to work:

e=eps(u)
sigm=sigma(e)
plot(sigm[0,0],mode='color')