Problem implementing effective stress function method with FEniCS

Hi all! I have been working on implementing the effective stress function method (Kojic and Bathe 1987) for plasticity-creep problem within the framework of FEniCS.

I followed the Numerical tours of continuum mechanics using FEniCS example (Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation) and studied the MFront implementation for the same problem.

I’m quite new to both FEM and FEniCS, and I believe I encountered a problem getting an working tangent stiffness matrix to work with the form “a_Newton_ESF = inner(eps(v), dot(Ct, eps(u_)))*dx”, where Ct is tangent stiffness matrix. With this implementation, the convergence is suboptimal, also the results seem wrong (nonsymmetric displacement although load should be symmetric along centerline of domain.)

In the Kojic and Bathe paper, Ct, to my understanding, is a 3 by 3 matrix. I used the above form. But I get very slow convergence and wrong results, especially with yielding. How does this 3 by 3 matrix (order-2 tensor) dot a strain (epsilon, another 3 by 3 matrix) give me stress? It seems that the tangent stiffness should be a order-4 tensor, right? Do you see what I do wrong here?

The code:

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=2, threshold = 1000,suppress=True,linewidth=220)
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)
sig0 = Constant(250.)  # yield strength
Et = E/100.  # tangent modulus
H = E*Et/(E-Et)  # hardening modulus
Ep = E*Et/(E-Et)  # hardening modulus, E_p in Anna_Pandolfi
aE = (1.0+nu)/E

Length = 1.0
Height = 1.0
Nx = 1
Ny = 1
mesh = RectangleMesh(Point(0,0), Point(Length, Height), Nx, Ny, 'left')

# Function spaces will involve a standard CG space for the displacement whereas internal
# state variables such as plastic strains will be represented using a **Quadrature** element.
# This choice will make it possible to express the complex non-linear material constitutive
# equation at the Gauss points only, without involving any interpolation of non-linear
# expressions throughout the element. It will ensure an optimal convergence rate
# for the Newton-Raphson method. See Chapter 26 of the `FEniCS book <>`_.
# We will need Quadrature elements for 4-dimensional vectors and scalars, the number
# of Gauss points will be determined by the required degree of the Quadrature element
# (e.g. ``degree=1`` yields only 1 Gauss point, ``degree=2`` yields 3 Gauss points and
# ``degree=3`` yields 6 Gauss points (note that this is suboptimal))::

deg_u = 2
deg_stress = 2

V = VectorFunctionSpace(mesh, "CG", deg_u)
Vstr = VectorFunctionSpace(mesh, "DG", 1, dim=4)

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)

Te = TensorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, shape = (3,3), quad_scheme='default')
T = FunctionSpace(mesh, Te)
T_ = TensorFunctionSpace(mesh, "DG", 1, shape = (3,3))

# .. note::
#  In older versions, it was possible to define Quadrature function spaces directly
#  using ``FunctionSpace(mesh, "Quadrature", 1)``. This is no longer the case since
#  FEniCS 2016.1 (see `this issue <>`_). Instead, Quadrature elements must first be defined
#  by specifying the associated degree and quadrature scheme before defining the
#  associated FunctionSpace.
# Various functions are defined to keep track of the current internal state and
# currently computed increments::
# Created functions in various function spaces.

sig = Function(W)
sig_old = Function(W)
u = Function(V, name="Total displacement")
u1 = Function(V, name="Total displacement")
tmp_u = Function(V, name="Total displacement")
tmp_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)
new_epl = Function(W)
neps = Function(W)

# Boundary conditions:

bctol = 1E-10
def left(x):
    return near(x[0],0., bctol)
def right(x):
    return near(x[0],Length, bctol)
def bottom(x):
    return near(x[1],0., bctol)
def top(x):
    return near(x[1],Height, bctol)

bc = [DirichletBC(V.sub(0), Constant(0.0), left), DirichletBC(V.sub(0), Constant(0.0), right),
        DirichletBC(V,Constant((0.0,0.0)), top)]

n = FacetNormal(mesh)
# Return symbolic facet normal for given mesh.
q_lim = float(2/sqrt(3)*ln(8.1)*sig0)
loading = Expression("q*t", q=q_lim, t=0, degree=2)

def F_ext(v):
    return loading*dot(n, v)*ds
# Variational form of the external load

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 v_eps(v):
    e = sym(grad(v))
    return as_vector([e[0, 0], e[1, 1], 0., e[0, 1]])

def sigma(eps_el):
    return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el

def as_3D_tensor(X):
    # convert a 4-d vector to 3d tensor representation
    return as_tensor([[X[0], X[3], 0], [X[3], X[1], 0], [0, 0, X[2]]])

def proj_sig_ESF(new_eps_vec, old_epl):
    # inputs are two vectors: new total strain and old plastic strain 
    new_eps_tensor = as_3D_tensor(new_eps_vec)

    dev_new_eps_tensor = dev(new_eps_tensor)
    mean_new_eps_tensor = tr(new_eps_tensor)/3.0
    epp = dev_new_eps_tensor - as_3D_tensor(old_epl) # calculate current epp

    # ESF parameter d is calculated with epp
    d = sqrt(1.5*inner(epp,epp))

    # Calculate trial equivalent stress
    sig_eq_trial = d/aE

    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)

    # compute new values of sigma (new_sig), delta plastic strain (delta_epl)
    new_sig_dev = 1.0/(aE+dlambda)*epp # new deviatoric stress (big S)
    new_sig_mean = E/(1.0-2*nu)*mean_new_eps_tensor
    new_sig = new_sig_dev + new_sig_mean*Identity(3)

    delta_epl = dlambda * new_sig_dev # increment of equivalent plastic strain

    new_epl = as_3D_tensor(old_epl) + delta_epl

    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_epl[0, 0], new_epl[1, 1], new_epl[2, 2], new_epl[0, 1]]), as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]])

# --------------------------------------------
# Global problem and Newton-Raphson procedure
# --------------------------------------------
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

a_Newton_ESF = inner(eps(v), dot(as_3D_tensor(Ct_vec), eps(u_)))*dxm
res_ESF = F_ext(u_) - inner(eps(u_), as_3D_tensor(sig))*dxm

# 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)
        # Solve local (cell-wise) problems A_e x_e = b_e where A_e and b_e are the cell element tensors.
        return u

# We now define the global Newton-Raphson loop. We will discretize the applied

Nitermax, tol = 30, 1e-13  # parameters of the Newton-Raphson procedure
Nincr = 1
load_steps = np.linspace(0, 1.0, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))

for (i, t) in enumerate(load_steps):
    loading.t = t
    nRes0 = 1.0
    nRes = 2.0
    print("Increment:", str(i+1))
    niter = 0
    while nRes > tol and niter < Nitermax: # NR loop starts

        local_project(v_eps(u1), W, neps)
        Ct_vec_, new_epl_, new_sig_ = proj_sig_ESF(neps, old_epl)
        #neps is vector, old_epl is vector
        local_project(new_epl_, W, new_epl)
        local_project(new_sig_, W, sig)
        local_project(Ct_vec_, W, Ct_vec)

        A, Res = assemble_system(a_Newton_ESF, res_ESF, bc)
        solve(A, du.vector(), Res, "mumps")


        A, Res = assemble_system(a_Newton_ESF, res_ESF, bc)
        nRes = Res.norm("l2")
        print("niter: ", niter," Residual:", nRes)
        niter += 1
    old_epl.assign(new_epl) # new

c = plot(sqrt(u[0]*u[0] + u[1]*u[1]))

Kojić, Miloš, and Klaus‐Jürgen Bathe. “The ‘effective‐stress‐function’algorithm for thermo‐elasto‐plasticity and creep.” International journal for numerical methods in engineering 24.8 (1987): 1509-1532.

I think maybe the confusion here is between Voigt notation and standard tensor notation. When Kojić and Bathe give a 3\times 3 material stiffness matrix, they are assuming it multiplies a vector of the unique plane-strain components,

\left\lbrack\begin{array}{c}\epsilon_{11}\\\epsilon_{22}\\2\epsilon_{12}\end{array}\right\rbrack\text{ ,}

where the factor of 2 on the shear component is because Voigt notation strains traditionally use the engineering shear strains, e.g., \gamma_{12} := 2\epsilon_{12}, (which represent changes in angle). The result of the matrix–vector product is then a vector of unique in-plane stress components,


(where there is no factor of two on the shear stress). You can use the as_vector and as_matrix functions in UFL to convert between Voigt and tensor notations, if you want to use formulas for material stiffness tensors given in Voigt notation.


Hi Prof. Kamensky,

Your reply is very helpful. Thank you very much! I updated the code that gives me the stress from

dot(as_3D_tensor(Ct_vec), eps(u_))


a_Newton_ESF = inner(eps(v), dot(as_3D_tensor(Ct_vec), eps(u_)))*dxm

to a function:

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., nu*(tmp_sig[0]+tmp_sig[1])]])
                                #I still want to build the stress tensor so that 
                                #I can do the same inner(tensor_1, tensor_2) 
                                #in a_Newton_ESF
    return tmp_sig_tensor


a_Newton_ESF = inner(eps(v), sigma_tang(as_3D_tensor(Ct_vec), eps(u_)))*dxm

Now I get results that converge faster (but not as fast as quadratic convergence, which is supposed to be for a properly implemented Newton method) and symmetric. But there is a roughly 1% relative difference between my displacement and the results from the Numerical tours of continuum mechanics using FEniCS example. I suspect that I still did something wrong.

One candidate is the sig_zz in the sigma tensor that I built. With the special form of the Kojić and Bathe tangent stiffness matrix,

is it still correct to write the following?


An example of the Convergence sequence look like this:

Increment: 20
niter:  0  Residual: 10.811880553254007
niter:  1  Residual: 4.083067613092935
niter:  2  Residual: 1.5445647408359051
niter:  3  Residual: 0.5922397032102841
niter:  4  Residual: 0.22815264530316065
niter:  5  Residual: 0.08804580745109687
niter:  6  Residual: 0.03400011136316638
niter:  7  Residual: 0.013132971097308292
niter:  8  Residual: 0.005073275508134279
niter:  9  Residual: 0.0019598843076917777
niter:  10  Residual: 0.0007571445613006976
niter:  11  Residual: 0.00029250254034672297
niter:  12  Residual: 0.00011300078003228369
niter:  13  Residual: 4.365495823447619e-05
niter:  14  Residual: 1.6864982172362073e-05
niter:  15  Residual: 6.515357271215003e-06
niter:  16  Residual: 2.5170448404153996e-06
niter:  17  Residual: 9.723944278945909e-07
niter:  18  Residual: 3.756605494327771e-07
niter:  19  Residual: 1.4512803216654037e-07
niter:  20  Residual: 5.606554481910209e-08
niter:  21  Residual: 2.1659816526934866e-08
niter:  22  Residual: 8.366728793717869e-09
niter:  23  Residual: 3.234467509830643e-09
niter:  24  Residual: 1.2495189461830964e-09
niter:  25  Residual: 4.815444576443768e-10
niter:  26  Residual: 1.8614615210521758e-10
niter:  27  Residual: 7.185429165405678e-11
niter:  28  Residual: 2.6934158761182686e-11
niter:  29  Residual: 1.0458476845403674e-11
niter:  30  Residual: 4.3039617281662785e-12
niter:  31  Residual: 2.7202713881723804e-12
niter:  32  Residual: 1.593918329340801e-12
niter:  33  Residual: 2.1834072240720706e-13

 u on nodes:
0 1.000000 1.000000 -0.000000 0.000000
1 1.000000 0.000000 -0.000000 -0.007495
2 0.000000 1.000000 0.000000 -0.000000
3 0.500000 0.500000 -0.000000 -0.003747
4 0.500000 1.000000 -0.000000 -0.000000
5 1.000000 0.500000 -0.000000 -0.003747
6 0.000000 0.000000 0.000000 -0.007495
7 0.000000 0.500000 0.000000 -0.003747
8 0.500000 0.000000 0.000000 -0.007495

It looks like this component of sigma_tang just gets multiplied by zero in the contraction with eps(v), so, if I’m tracing the code correctly, it shouldn’t influence the result.

That makes sense! I also tried to arbitrarily modify the sigma_zz term and indeed there is no effect.