Hi, I am trying to build a displacement-controlled version of this tutorial Elasto-plastic analysis of a 2D von Mises material in which the author has used external force.
I tried to follow these posts Plasticity in fenics and Custom Newton Solver problem with Dirichlet conditions where @bleyerj states that to make a displacement controlled model, we need to keep 2 things in mind -
- Add a DirichletBC with imposed value equal to the displacement increment over the time step and not the total displacement.
- Only the first iteration should have a non-zero DirichletBC all other iterations must be computed with zero DirichletBC, otherwise you will accumulate imposed displacements over all iterations and will not converge.
Following both advices, I modified the loop like this -
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, L, bc_u)
for bc in bc_u:
bc.apply(du.vector())
bc.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0
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, L, bc_u)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)
But the solver diverges, I think I’ve not understood both the points due to which there is divergence. Therefore, I request to kindly give me suggestions on the same.
Thanks.
PS: Here is my MWE-
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
# elastic parameters
E = Constant(70e3 )
nu = 0.3
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
Et = E/100. # tangent modulus
H = E*Et/(E-Et) # hardening modulus
sig0 = Constant(250.)
from mshr import *
domain= Rectangle(Point(0., 0.), Point(1., 1.))
mesh = generate_mesh(domain, 20)
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)
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
# Boundary conditions
class top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1., 1e-8)
class bot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)
boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
top().mark(boundaries, 1)
bot().mark(boundaries, 2)
loading = Expression("t", t = 0.0, degree=1)
bc_top = DirichletBC(V.sub(1), loading, boundaries, 1)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)
bc_u = [bc_bottom, bc_top]
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 = inner(eps(v), sigma_tang(eps(u_)))*dxm
L = -inner(eps(u_), as_3D_tensor(sig))*dxm
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
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, L, bc_u)
for bc in bc_u:
bc.apply(du.vector())
bc.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0
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, L, bc_u)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)