Hello, following this 2D elastoplastic example, I am trying to modify the code to a 3D problem, I believe I have set the dimensions correctly to full 3D but the Newton Raphson doesn’t converge. in the following code I try to point out things that I had changed (other than domain and BCs, and loading force), while the rest follows the tutorial. I have tried to do a test by setting the yield stress sigma0 to a very large number( = dl.Constant(1e9.) for example), to keep the simulation in the elastic region and still having the same problem of Newton Raphson not converging, can anyone help spotting any issue in the code?
(update: the same mesh and BCs and load were testing in a separate hyperplasticity code and were fine, so I don’t think the problem would be in there)
##### domain and BCs stuff
import dolfin as dl
import numpy as np
dl.parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
import matplotlib.pyplot as plt
longdim=1000.
shortdim=150.
mesh = dl.BoxMesh(dl.Point(-longdim/2., shortdim/2., -shortdim/2.), dl.Point(longdim/2., -shortdim/2., shortdim/2.), 10, 4, 4)
class LeftSurf(dl.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and dl.near(x[0], -longdim/2.)
class RightSurf(dl.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and dl.near(x[0], longdim/2.)
boundary_parts = dl.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary_parts.set_all(0) #marks whole cube as domain 0
left = LeftSurf()
left.mark(boundary_parts, 1)
right = RightSurf()
right.mark(boundary_parts, 2)
ds=dl.Measure('ds', domain=mesh,subdomain_data=boundary_parts)
### following tutorial
E = dl.Constant(70e3)
nu = dl.Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = dl.Constant(250.)
Et = E/100.
H = E*Et/(E-Et)
deg_u = 2
deg_stress = 2
V = dl.VectorFunctionSpace(mesh, "CG", deg_u)
u = dl.Function(V, name="Total displacement")
du = dl.Function(V, name="Iteration correction")
Du = dl.Function(V, name="Current increment")
v = dl.TrialFunction(V)
u_ = dl.TestFunction(V)
##### here the dimension of We is updated to 6 corresponding
#####to the 6 unique components in Voigt notation in 3D
We = dl.VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=6, quad_scheme='default')
W = dl.FunctionSpace(mesh, We)
### following tutorial
sig = dl.Function(W)
sig_old = dl.Function(W)
n_elas = dl.Function(W)
W0e = dl.FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = dl.FunctionSpace(mesh, W0e)
beta = dl.Function(W0)
p = dl.Function(W0, name="Cumulative plastic strain")
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dl.dx(metadata=metadata)
##### some more BCs and loading stuff
bc = [dl.DirichletBC(V, dl.Constant((0.,0.,0.)), LeftSurf()) , dl.DirichletBC(V, dl.Constant((0.,0.,0.)), RightSurf())]
class PointLoad_t(dl.UserExpression):
def __init__(self,vl, t, **kwargs):
super().__init__(**kwargs)
self.value = vl
self.tol = 1.
self.t = t
def eval(self, values, x):
if dl.near(x[0],0.,self.tol) and dl.near(x[1],0.,self.tol) and dl.near(x[2],-shortdim/2.,self.tol):
values[0] = 0.
values[1] = 0.
values[2] = self.value * self.t
else:
values[0] = 0.
values[1] = 0.
values[2] = 0.
def value_shape(self):
return (3,)
loading = PointLoad_t( 1., t=0., degree=1 )
def F_ext(v):
return dl.dot(loading, v)*dxm
#### updated vector to tensor conversion with 6 Voigt components
def as_3D_tensor(X):
return dl.as_tensor([[X[0], X[3], X[4]],
[X[3], X[1], X[5]],
[X[4], X[5], X[2]]])
def eps(v):return dl.sym(dl.grad(v))
def sigma(eps_el):
return lmbda*dl.tr(eps_el)*dl.Identity(3) + 2*mu*eps_el
ppos = lambda x: (x+abs(x))/2.
#### output Voigt vectors are updated to 6 components
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dl.dev(sig_elas)
sig_eq = dl.sqrt(3/2.*dl.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 dl.as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1], new_sig[0, 2], new_sig[1, 2]]), \
dl.as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1], n_elas[0, 2], n_elas[1, 2]]), \
beta, dp
#### the rest of the code follows tutorial
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*dl.inner(N_elas, e)*N_elas-2*mu*beta*dl.dev(e)
def local_project(v, V, u=None):
dv = dl.TrialFunction(V)
v_ = dl.TestFunction(V)
a_proj = dl.inner(dv, v_)*dxm
b_proj = dl.inner(v, v_)*dxm
solver = dl.LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = dl.Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
file_results = dl.XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = dl.FunctionSpace(mesh, "DG", 0)
p_avg = dl.Function(P0, name="Plastic strain")
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))
a_Newton = dl.inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -dl.inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)
for (i, t) in enumerate(load_steps):
loading.t = t
A, Res = dl.assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(dl.Constant((0., 0., 0.))). ##### updated to 3 components
print("Increment:", str(i+1), ', t = ', t)
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
dl.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 = dl.assemble_system(a_Newton, res, bc)
nRes = Res.norm("l2")
print("Iter #", str(niter) ," Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
file_results.write(u, t)
p_avg.assign(dl.project(p, P0))
file_results.write(p_avg, t)
results[i+1, :] = (u(0.,0., -shortdim/2.)[3], t)