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')
plt.ion()
plot(mesh)
plt.show()
plt.pause(1.5)
plt.ioff()
# 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 <https://fenicsproject.org/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 <https://bitbucket.org/fenics-project/dolfin/issues/757/functionspace-mesh-quadrature-1-broken-in>`_). 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)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
# Solve local (cell-wise) problems A_e x_e = b_e where A_e and b_e are the cell element tensors.
return u
else:
solver.solve_local_rhs(u)
return
# 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
u1.assign(u)
print("Increment:", str(i+1))
niter = 0
while nRes > tol and niter < Nitermax: # NR loop starts
u1.assign(u1+du)
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")
u.assign(u1)
old_epl.assign(new_epl)
A, Res = assemble_system(a_Newton_ESF, res_ESF, bc)
nRes = Res.norm("l2")
print("niter: ", niter," Residual:", nRes)
niter += 1
u.assign(u+Du)
sig_old.assign(sig)
old_epl.assign(new_epl) # new
plt.ion()
c = plot(sqrt(u[0]*u[0] + u[1]*u[1]))
plt.colorbar(c)
plt.show()
plt.pause(3)
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.
