Hi, sorry for the late reply, I tried to modify your code according to best of my knowledge to convert it into displacement controlled as stated in this post above and in the post Plasticity in fenics both time by @bleyerj but faced the same issue which is NewtonSolver converges first then at 3rd or 4th iteration gives nan value.
The modified code -
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)
# material 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
H = E*Et/(E-Et) # hardening modulus
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)
# 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)
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(V, Constant((0.0,0.0)), boundaries, 2)
bctop = DirichletBC(V.sub(1), load, boundaries, 1)
bc = [bcbot, bctop]
n = FacetNormal(mesh)
ds = Measure("ds")(subdomain_data=boundaries)
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)
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_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -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
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
u_r = 0.05
Nitermax, tol = 5, 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):
load.t = t*u_r
A, Res = assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "lu")
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_Newton, res, bc)
for bc_i in bc:
bc_i.homogenize()
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)
Gives the output -
Increment: 1
Residual: 8193.78486479144
Residual: 12.584295265311695
Residual: 27.48110701737934
Residual: 2.258960298045684
Residual: 0.7402582538956958
Increment: 2
Residual: 0.8915496012667151
Residual: 0.44893033099816865
Residual: 0.07483532946744006
Residual: 0.009844688153043523
Residual: 0.0007199674903077572
Increment: 3
Residual: 0.0002757595008993816
Residual: 2.2939673222581952e-05
Residual: 8.785003486230311e-07
Residual: 6.257779286847474e-09
Residual: 8.034401067089978e-14
Increment: 4
Residual: 7.378512786780972e-14
Residual: nan
Increment: 5
Increment: 6
Increment: 7
Increment: 8
Increment: 9
Increment: 10
Increment: 11
Increment: 12
Increment: 13
Increment: 14
Increment: 15
Increment: 16
Increment: 17
Increment: 18
Increment: 19
Increment: 20
And one more thing I am not sure what you were trying to do in these lines of code -
ST = TensorFunctionSpace(mesh, "DG", 0, shape =(3,3))
sigma_t = Function(P0)
load_vertical = Function(P0, name="vertical reaction")
sigma_tensor = as_3D_tensor(sig_)
sigma_t.assign(local_project(sig_[1], P0))
normal = as_vector((n[0],n[1], 0))
Traction = dot(sigma_t, normal)
f_y = Traction[1] *ds(1)
file_results.write(str(assemble(f_y)), t)
Because running it yields the error -
Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
/tmp/ipykernel_170/2372490712.py in <module>
41 sigma_t.assign(local_project(sig_[1], P0))
42 normal = as_vector((n[0],n[1], 0))
---> 43 Traction = dot(sigma_t, normal)
44 f_y = Traction[1] *ds(1)
45 file_results.write(str(assemble(f_y)), t)
/usr/lib/python3/dist-packages/ufl/operators.py in dot(a, b)
174 if a.ufl_shape == () and b.ufl_shape == ():
175 return a * b
--> 176 return Dot(a, b)
177
178
/usr/lib/python3/dist-packages/ufl/tensoralgebra.py in __new__(cls, a, b)
188 # Checks
189 if not ((ar >= 1 and br >= 1) or scalar):
--> 190 error("Dot product requires non-scalar arguments, "
191 "got arguments with ranks %d and %d." % (ar, br))
192 if not (scalar or ash[-1] == bsh[0]):
/usr/lib/python3/dist-packages/ufl/log.py in error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))
159
160 def begin(self, *message):
UFLException: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.
which makes sense because sigma_t is scalar.