Contact Mechanics - 2D problem

variational formulation
1

Hi,

I am trying to solve a 2D contact problem based on the 3D code available here:

But I can’t figure out why my code is not running. My code is the following

Ny=30 # Refinamento da malha na direção y

Nx=30 # Refinamento da malha na direção x

mesh = fn.RectangleMesh(fn.Point(0, -1), fn.Point(1, 0),Nx,Ny) # gera a malha 

mesh.coordinates()[:,0]=mesh.coordinates()[:,0]**2

mesh.coordinates()[:,1]=-mesh.coordinates()[:,1]**2

fn.plot(mesh, title="Malha")

class Top(fn.SubDomain):

    def inside(self, x, on_boundary):

        return fn.near(x[1], 0.) and on_boundary

def symmetry_x(x, on_boundary):

        return fn.near(x[0], 0.) and on_boundary

def bottom(x, on_boundary):

        return fn.near(x[1], -1.) and on_boundary

# exterior facets MeshFunction

facets = fn.MeshFunction("size_t", mesh, 1)

facets.set_all(0)

Top().mark(facets, 1)

ds = fn.Measure('ds', subdomain_data=facets)

R = 0.5

d = 0.02

obstacle = fn.Expression("-d+pow(x[0],2)/2/R", d=d, R=R, degree=2)

V = fn.VectorFunctionSpace(mesh, "CG", 1)

V2 = fn.FunctionSpace(mesh, "CG", 1)

V0 = fn.FunctionSpace(mesh, "DG", 0)

u = fn.Function(V, name="Displacement")

du = fn.TrialFunction(V)

u_ = fn.TestFunction(V)

bc =[fn.DirichletBC(V, fn.Constant((0., 0.)), bottom),

     fn.DirichletBC(V.sub(0), fn.Constant(0.), symmetry_x)]

E = fn.Constant(10.)

nu = fn.Constant(0.3)

mu = E/2/(1+nu)

lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

    return fn.sym(fn.grad(v))

def sigma(v):

    return lmbda*fn.tr(eps(v))*fn.Identity(2) + 2.0*mu*eps(v)

def ppos(x):

    return (x+np.abs(x))/2.

pen = fn.Constant(1e4)

form = fn.inner(sigma(u), eps(u_))*fn.dx + pen*fn.dot(u_[1], ppos(u[1]-obstacle))*fn.ds(1)

J = fn.derivative(form, u, du)

problem = fn.NonlinearVariationalProblem(form, u, bc, J=J)

solver = fn.NonlinearVariationalSolver(problem)

solver.parameters["newton_solver"]["linear_solver"] = "cg"

solver.parameters["newton_solver"]["preconditioner"] = "ilu"

solver.solve()

import matplotlib.pyplot as plt

p = fn.plot(sigma(u)[1,1], mode='color') 

plt.colorbar(p)

plt.title(r"$\sigma_{yy}$",fontsize=26)

plt.show()

I tried different nonlinear solvers but the result is the same. What is wrong?

2

Could you please supply the error message you are getting, as well as specifying what version of fenics you are using.

3

The difference comes from the fact that the tutorial is using a wildcard import (i.e., from dolfin import *). The problem is that you set subdomain_data in ds, but then use fn.ds in form. You can either make the change

fn.ds = fn.Measure('ds', subdomain_data=facets)

(where I’ve added an fn. at the beginning) or use ds(1) (without the fn.) to integrate the boundary term.

2 Likes
4

That’s true. Thank you so much. It solved my problem. Now it works.

Cheers.

5

Hey, how would one model the indenter explicitly instead of using the Expression without any mesh?

6

I am interested in this problem as well. You need to discretize both bodies and consider some penalty/augmented Lagrange formulation in order to impose displacement constraints and force equilibrium on the contact surfaces of both bodies. I don’t believe that there is any example available in the web and I don’t know if it is possible to solve this problem with Fenics yet. The best discussion that I found in the web was this:

https://fenicsproject.org/qa/880/surface-contact/

It was from many years ago (2013). At that time, to solve this problem was very difficult or even impossible with Fenics. As far as I know, few things have changed since then.

7

I discretized both bodies and also marked them differently as done in the fenics tutorial. My idea was to define a distance function between the lower boundary of the upper geometry and the upper boundary of the bottom geometry and use it for the penalty approach. But I am not sure how to do that. My problem is that I do not understand how to use the boundary values of the two boundaries as a function input.

8

I suggest you to have a look at this link:

https://dealii.org/current/doxygen/deal.II/step_42.html

It is not Fenics but it is well documented and has a code with the implementation. It is the best finite element code applied to contact mechanics that I found. I am trying to follow a similar approach in order to implement it in Fenics. But I am still in the beginning.

9

Hi,
i want to modify the points of mesh as
mesh = fn.RectangleMesh(fn.Point(0, 0), fn.Point(1, 1),Nx,Ny) # gera a malha ,
so that all the coordinates are >= 0. For that i still have to make 4 modification of the code as

import fenics as fn

Ny=30 # Refinamento da malha na direção y

Nx=30 # Refinamento da malha na direção x

###################  Modification1  #############
#mesh = fn.RectangleMesh(fn.Point(0, -1), fn.Point(1, 0),Nx,Ny) # gera a malha 
mesh = fn.RectangleMesh(fn.Point(0, 0), fn.Point(1, 1),Nx,Ny) # gera a malha 
###################  Modification1  #############


mesh.coordinates()[:,0]=mesh.coordinates()[:,0]**2


###################  Modification2  #############
#mesh.coordinates()[:,1]=-mesh.coordinates()[:,1]**2
mesh.coordinates()[:,1]=mesh.coordinates()[:,1]**2
###################  Modification2  #############


#fn.plot(mesh, title="Malha")

class Top(fn.SubDomain):

    def inside(self, x, on_boundary):

###################  Modification3  #############        
 #        return fn.near(x[1], 0.) and on_boundary
        return fn.near(x[1], 1.) and on_boundary
###################  Modification3  #############


def symmetry_x(x, on_boundary):

        return fn.near(x[0], 0.) and on_boundary

def bottom(x, on_boundary):

###################  Modification4  #############
 #        return fn.near(x[1], -1.) and on_boundary
        return fn.near(x[1], 0.) and on_boundary
###################  Modification4  #############



# exterior facets MeshFunction

facets = fn.MeshFunction("size_t", mesh, 1)

facets.set_all(0)

Top().mark(facets, 1)

fn.ds = fn.Measure('ds', subdomain_data=facets)

R = 0.5

d = 0.02

###################  Modification5  #############
#obstacle = fn.Expression("-d+pow(x[0],2)/2/R", d=d, R=R, degree=2)
obstacle = fn.Expression("1.0-d+pow(x[0],2)/2/R", d=d, R=R, degree=2)
###################  Modification5  #############


V = fn.VectorFunctionSpace(mesh, "CG", 1)

V2 = fn.FunctionSpace(mesh, "CG", 1)

V0 = fn.FunctionSpace(mesh, "DG", 0)

u = fn.Function(V, name="Displacement")

du = fn.TrialFunction(V)

u_ = fn.TestFunction(V)

bc =[fn.DirichletBC(V, fn.Constant((0., 0.)), bottom),

     fn.DirichletBC(V.sub(0), fn.Constant(0.), symmetry_x)]

E = fn.Constant(10.)

nu = fn.Constant(0.3)

mu = E/2/(1+nu)

lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

    return fn.sym(fn.grad(v))

def sigma(v):

    return lmbda*fn.tr(eps(v))*fn.Identity(2) + 2.0*mu*eps(v)

def ppos(x):

    return (x+np.abs(x))/2.

pen = fn.Constant(1e4)

form = fn.inner(sigma(u), eps(u_))*fn.dx + pen*fn.dot(u_[1], ppos(u[1]-obstacle))*fn.ds(1)

J = fn.derivative(form, u, du)

problem = fn.NonlinearVariationalProblem(form, u, bc, J=J)

solver = fn.NonlinearVariationalSolver(problem)

solver.parameters["newton_solver"]["linear_solver"] = "cg"

solver.parameters["newton_solver"]["preconditioner"] = "ilu"

solver.solve()

import matplotlib.pyplot as plt

p = fn.plot(u, mode="displacement") 

plt.colorbar(p)

plt.show()

But the output is that all displacements are 0. It is very strange^^

10

Sorry for the delay in answering your question. What is wrong is Modification 5. Just return to the original one and everything will work. This means the imposed displacement. So, it should be near zero.

11

Thanks a lot for your reply.
Because in your example the mesh points are “fn.Point(0, -1), fn.Point(1, 0)”, but in my example are “fn.Point(0, 0), fn.Point(1, 1)”, so the top of the mesh in your example is y=0 and in my example is y=1, therefore in Modification 5 the obstacle should be added by 1.0.

12

No, obstacle can be seen as the penetration of the rigid body into the elastic body. It should be a small number (absolute value), much less than 1. The largest absolute value of penetration is at x[0]=0, where it is equal to -d. If the obstacle is positive, it means no contact.

14

Oh you are right. Thanks a lot ^^