Implementing effective stress function method for plasticity

Hi all! I have been working on implementing the effective stress function method (Kojić and Bathe 1987) for plasticity-creep problem within the framework of FEniCS. Previously I have posted a question regarding this topic and have got very helpful feedback. After a while I am back to this problem. I have difficulties achieving the quadratic convergence. I wish I can get some ideas on what I might do wrong.

I followed the Numerical tours of continuum mechanics using the famous FEniCS example by Jeremy Bleyer (Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation).

The MWE is attached, note to run it, you will need to download the mesh file from the above online example: mesh files and original code by Jeremy Bleyer.

The MWE simplified the original example to 2 load steps, 1st one elastic, and 2nd one would cause plasticity. However for the plastic load, the convergence is slow and is nowhere near the quadratic convergence, which suggests that I might get the tangent stiffness matrix wrong.

The tangent stiffness matrix is updated in the update_sig function, which takes the current displacement and the plastic strain from last step as input. From Kojić and Bathe, in the case with only plasticity (and bilinear stress-strain relationship), the effective stress can be solved without bisection as

where d is defined as


e’ is the deviatoric strain.

Again with bilinear stress-strain relationship, the plasticity-related parameter depends only on the effective stress

The computational procedure starts by calculating the elastic solution

If yielding happened, the effective stress is taken as eq. (44).

With the effective stress solved, delta lambda can be solved, and the tangent stiffness matrix is updated as

Can you see what I did wrong in the code? (I used FEniCS conditional which could be tricky…) Or did I misunderstand the algorithm?

from dolfin import *
import numpy as np
import math
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
Ep = E*Et/(E-Et)  # hardening modulus
aE = (1.0+nu)/E

Re, Ri = 1.3, 1.   # external/internal radius
mesh = Mesh("thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml")
ds = Measure('ds')[facets]

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)

# Various functions are defined to keep track of the current internal state and
# currently computed increments::

sig = Function(W)
sig_old = Function(W)
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)

Ct_vec = Function(W)
old_epl = Function(W)

# Boundary conditions

bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]

n = FacetNormal(mesh)
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)
loading = Expression("-q*t", q=q_lim, t=0, degree=2)

def F_ext(v):
    return loading*dot(n, v)*ds(4)

# -----------------------------
# Constitutive relation update
# -----------------------------

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):
    # convert dimension 4 vector to 3D tensor
    return as_tensor([[X[0], X[3], 0],
                      [X[3], X[1], 0],
                      [0, 0, X[2]]])

def update_sig(u, old_epl):
    new_eps = eps(u) # 3D tensor
    new_eps_mean = tr(new_eps)/3.0 * Identity(3) # 3D tensor
    new_eps_dev = new_eps - new_eps_mean # 3D tensor
    new_epp = new_eps_dev - as_3D_tensor(old_epl) # 3D tensor

    d = sqrt(1.5*inner(new_epp, new_epp))
    sig_eq_trial = d/(aE)

    # effective stress can be calculated without bisection for bilinear case
    new_sig_eq = conditional(lt(sig_eq_trial, sig0), sig_eq_trial, (2.0*Ep*d+3.0*sig0)/(2.0*Ep*aE+3.0))
    dlambda = conditional(lt(sig_eq_trial, sig0), Constant(0.), (3.0/2.0)*(1.0-sig0/new_sig_eq)/Ep)

    # next calculate new deviatoric stress from the new effective stress
    new_sig_dev = 1.0/(aE+dlambda)*new_epp
    new_sig_mean = E/(1.0-2*nu)*new_eps_mean
    new_sig = new_sig_dev + new_sig_mean

    delta_epl = dlambda * new_sig_dev # increment of equivalent plastic strain
    new_epl = as_3D_tensor(old_epl) + delta_epl

    new_sig = new_sig_dev + new_sig_mean

    Cp = 1.0/(aE + dlambda)
    Cm = E/(1.0-2.0*nu)
    C11 = 1.0/3.0*(2.0*Cp+Cm)
    C22 = C11
    C12 = 1.0/3.0*(Cm-Cp)
    C33 = 0.5*Cp

    return as_vector([C11, C22, C33, C12]), as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
           as_vector([new_epl[0, 0], new_epl[1, 1], new_epl[2, 2], new_epl[0, 1]])

def sigma_tang(Ct, e):
    tmp_e = as_vector([e[0, 0], e[1, 1], 2.0*e[0, 1]]) #Voigt notation of strain
    tmp_sig = dot(Ct, tmp_e) #Calculated stress in Voigt notation
    tmp_sig_tensor = as_tensor([[tmp_sig[0], tmp_sig[2], 0.], \
                                [tmp_sig[2], tmp_sig[1], 0.], \
                                [0., 0., 0.]]) #stress tensor
    return tmp_sig_tensor

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

metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

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

# The consitutive update defined earlier will perform nonlinear operations on
# the stress and strain tensors. These nonlinear expressions must then be projected
# back onto the associated Quadrature spaces. Since these fields are defined locally
# in each cell (in fact only at their associated Gauss point), this projection can
# be performed locally. For this reason, we define a ``local_project`` function
# that use the ``LocalSolver`` to gain in efficiency (see also :ref:`TipsTricksProjection`)
# for more details::

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)
    if u is None:
        u = Function(V)
        return u

# Before defining the Newton-Raphson loop, we set up the output file and appropriate
# FunctionSpace (here piecewise constant) and Function for output of the equivalent
# plastic strain since XDMF output does not handle Quadrature elements::

file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

# We now define the global Newton-Raphson loop.

Nitermax, tol = 10, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 2
load_steps = np.linspace(0, 0.9, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
    loading.t = t

    # Initial guess for tangent stiffness matrix
    Cp = 1.0/(aE)
    Cm = E/(1.0-2.0*nu)
    C11 = 1.0/3.0*(2.0*Cp+Cm)
    C22 = C11
    C12 = 1.0/3.0*(Cm-Cp)
    C33 = 0.5*Cp
    Ct_vec_ = as_vector([C11, C22, C33, C12])
    local_project(Ct_vec_, W, Ct_vec)

    A, Res = assemble_system(a_Newton, res, bc)
    nRes0 = Res.norm("l2")
    nRes = nRes0
    Du.interpolate(Constant((0, 0)))
    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")
        Ct_vec_, sig_, new_epl_ = update_sig(u+Du, old_epl) # input is current DeltaU (Du) and old plastic strain 
        local_project(Ct_vec_, W, Ct_vec)
        local_project(sig_, W, sig)
        A, Res = assemble_system(a_Newton, res, bc)
        nRes = Res.norm("l2")
        print("    Residual:", nRes)
        niter += 1
    old_epl.assign(local_project(new_epl_, W))

# ----------------
# Post-processing
# ----------------
    file_results.write(u, t)
    file_results.write(Du, t)