Hello everyone,
in the following MWE I just want to solve the displacement of a rectangle and is based on the following example: Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation
The code works fine when I impose a force on the top edge. However the same code does not converge when I try to define a displacement controlling dirichlet BC (here called “bctopy”) instead of the force.
MWE:
from dolfin import *
from mshr import *
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)
domain = Rectangle(Point(0,0), Point(1,1))
nxy = 40
mesh = generate_mesh(domain, nxy)
# elastic parameters
kappa = Constant(71245)
mu = Constant(27120)
lambda_= Constant(kappa - 2./3.*mu)
V = VectorFunctionSpace(mesh, "CG", 2)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=2, quad_scheme='default')
W = FunctionSpace(mesh, We)
sig = Function(W)
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
du_ = TrialFunction(V)
v = TestFunction(V)
# Boundary conditions for mesh D2
class bottom(SubDomain):
def inside(self,x,on_boundary):
tol = 1e-13
return abs(x[1]) < tol and on_boundary
class top(SubDomain):
def inside(self,x,on_boundary):
tol = 1e-13
return abs(x[1]-1) < tol and on_boundary
Top = top()
Bottom = bottom()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
Bottom.mark(boundaries,1)
Top.mark(boundaries,2)
loading = Expression("t", t = 0.0, degree=1)
bcbot = DirichletBC(V, Constant((0.0,0.0)), boundaries, 1)
bctopx = DirichletBC(V.sub(0), Constant(0.0), boundaries, 2)
bctopy = DirichletBC(V.sub(1), loading, boundaries, 2)
bc_u = [bcbot, bctopx, bctopy]
n = FacetNormal(mesh)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
def F_ext(v):
return loading*dot(n, v)*ds(2)
# -----------------------------
# Constitutive relation update
# -----------------------------
def eps(v):
return sym(grad(v))
def sigma(eps_el):
return lambda_*tr(eps_el)*Identity(2) + 2*mu*eps_el
sig = sigma(eps(u))
# --------------------------------------------
# Global problem and Newton-Raphson procedure
# --------------------------------------------
metadata = {"quadrature_degree": 2, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
a_Newton = inner(eps(v), sigma(eps(du_)))*dxm
res = -inner(eps(v), sig)*dxm# + F_ext(v)
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
Nitermax, tol = 50, 1e-8 # parameters of the Newton-Raphson procedure
Uincr = 9E-2
Nincr = 10
j = 0
t = 0
while j < Nincr:
j += 1
t += Uincr
loading.t = t
A, Res = assemble_system(a_Newton, res, bc_u)
nRes0 = Res.norm("l2")
nRes = nRes0
print("Increment:", str(j))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
u.assign(u+du)
A, Res = assemble_system(a_Newton, res, bc_u)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
file_results.write(u, t)
Dolfin version (installed via Docker): 2019.1.0
Ubunutu version: 20.04
Does anyone have an idea what the mistake is?
Your help is much appreciated.
Best,
Dejan