Dear all,
I am trying to implement elastoplasticity with contact boundary condition, followed from the examples provided here (elastoplasticity and contact). I have added the neuman boundary condition due to contact to the residue of the elastoplasticity weak form. Upon assembling, I always get a 0 residue, meaning no external forces on the domain. I believe the weak form is somehow wrongly defined, but I am not sure what the issue is.
The problem I am trying to solve is a single punch test with a cylindrical sample being punched from the top.
If it is not an issue with the weak form, what else could it be? Any help is greatly appreciated. I will attach the code below, and I am using fenics version 2019.2.0.dev0.
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
from mshr import *
def eps(v):
return sym(grad(v))
def sigma(eps_el):
return lmbda*tr(eps_el)*Identity(len(u)) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], X[4]],
[X[3], X[1], X[5]],
[X[4], X[5], 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], new_sig[0, 2], new_sig[1, 2]]), \
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
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)
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
class X_Plane(SubDomain):
def __init__(self, x1):
super().__init__()
self.x1 = x1
def inside(self, x, on_boundary):
return on_boundary and near(x[2], self.x1, tol)
# elastic parameters
E = Constant(2100e3)
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
tol = 1e-8
# Making Mesh (40 corresponds to the mesh density)
mesh = generate_mesh(Cylinder(Point(0, 0, 5), Point(0, 0, 0), 7.5, 7.5), 10)
# Save the mesh
File("Result/cylinder.pvd") << mesh
n = FacetNormal(mesh)
deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=6, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facets.set_all(0)
top = X_Plane(5)
top.mark(facets, 2)
ds = Measure('ds', subdomain_data=facets)
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)
bc_iter = [DirichletBC(V, Constant((0, 0, 0)), 'near(x[2], 0, 1e-8)')]
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
R = 1.25
d = 0.02
obstacle = Expression("-d+(pow(x[0],2)+pow(x[1], 2))/2/R", d=d, R=R, degree=2)
pen = Constant(1e4)
def F_ext(u_test, u_trial):
return pen*dot(u_test[2], ppos(u_trial[2] - obstacle))*ds(2)
a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
v = Function(V)
res = -inner(eps(u_), as_3D_tensor(sig))*dxm + pen*dot(u_[2], ppos(v[2] - obstacle))*ds(2)
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
disp = File('Result/displacement.pvd')
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))
u_inf_max = 0.1
delta_u = 1e-04
u.interpolate(Constant((0, 0, 0)))
p.interpolate(Constant(0))
step = 0
for obstacle.d in np.linspace(2e-2, 2e-1, num=100):
step += 1
print('===================================================')
print(' Total steps: ' + str(step))
print('===================================================')
obst = interpolate(obstacle, FunctionSpace(mesh, 'CG', 1))
A, Res = assemble_system(a_Newton, res, bc_iter)
print(np.max(Res[:]))
print(np.min(Res[:]))
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0, 0)))
niter = 0
# print(nRes)
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res)
u.assign(u+Du)
# File('displacement.pvd') << u
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_iter)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
disp << u
print('Simulation completed')