Hi
I am trying to implement a contact problem for 3d solid mechanics. I used this example ( 3D Solid Mechanics — NGS-Py 6.2.2405 documentation (ngsolve.org)) with few changes in the geometry, and added 3d contact problem from this example ( 6.2 Contact Problems — NGS-Py 6.2.2405 documentation (ngsolve.org)), but the solver did not converge. Can someone please tell me what I am doing wrong?
Here is my code:
from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.krylovspace import CGSolver
box1 = Box((0,0,0), (1,3,0.6))
box2 = Box((-1,0,0), (2,-1,0.6))
box1.faces.Min(Y).name = "contact1"
box2.faces.Max(Y).name = "contact2"
geo = Compound([box1, box2])
geo.faces.Min(Y).name = "fix"
geo.faces.Max(Y).name = "force"
mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1))
E, nu = 210, 0.2
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
def Stress(strain):
return 2*mu*strain + lam*Trace(strain)*Id(3)
fes = VectorH1(mesh, order=3, dirichlet="fix")
u,v = fes.TnT()
gfu = GridFunction(fes)
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
pre = Preconditioner(a, "bddc")
a.Assemble()
force = CF( (0,1e-3,0) )
f = LinearForm(force*v*ds("force")).Assemble()
contact = ContactBoundary(mesh.Boundaries("contact1"), mesh.Boundaries("contact2"))
X = CoefficientFunction((x,y,z))
cf = (X + u - (X.Other() + u.Other())) * specialcf.normal(3)
contact.AddIntegrator(IfPos(cf, 1e14*cf*(v-v.Other())*contact.normal, 0))
with TaskManager():
contact.Update(gfu, a)
inv = CGSolver(a.mat, pre, printrates=True, tol=1e-8)
gfu.vec.data = inv * f.vec
with TaskManager():
fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
gfstress = GridFunction(fesstress)
gfstress.Interpolate (Stress(Sym(Grad(gfu))))
print(gfstress(mesh(0,2,0))[4])
Draw (Norm(gfstress), mesh, deformation=4*gfu, draw_vol=False, order=3)