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
and
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)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
# 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")
Du.assign(Du+du)
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
u.assign(u+Du)
sig_old.assign(sig)
old_epl.assign(local_project(new_epl_, W))
# ----------------
# Post-processing
# ----------------
file_results.write(u, t)
file_results.write(Du, t)