Contact problem for 3d solid mechanics with CGSolver

NGSolve
1

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)
1 Like