Is there a contact model in FEniCS for two imported physical volumes?

[michi-ost]

Hi there
I am really new to fenics, but i worked the tutorials through. Now I am trying to do some more stuff to get into fenics.

https://1drv.ms/u/s!AocekCrLgWbAihmCF7NaIQBNEIit?e=HEuUJC

As you can see in the attached image, I try to simulate a pulling force on a physical surface (2) of a geometry. The second geometry is fixed at its left and right physical planes (1 + 3).
The code I created works but I do not know how to define a specific Force on the surface 2. In addition I have not found any contact body model. Can someone help me to get this to bodies connected?

I found this: https://comet-fenics.readthedocs.io/en/latest/demo/contact/penalty.html but I am not sure how to implement this in my code.
I will upload the .geo code in separat code.

Thank you very much.

import meshio

import matplotlib.pyplot as plt
import numpy as np
from fenics import *
from ufl import nabla_div
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp, Measure,
DirichletBC, FunctionSpace, Constant, TrialFunction,
TestFunction, dot, grad, dx, Function, solve

msh = meshio.read("./winkel.msh")

for cell in msh.cells:
if cell.type == “triangle”:
triangle_cells = cell.data
elif cell.type == “tetra”:
tetra_cells = cell.data

for key in msh.cell_data_dict[“gmsh:physical”].keys():
if key == “triangle”:
triangle_data = msh.cell_data_dict[“gmsh:physical”][key]
elif key == “tetra”:
tetra_data = msh.cell_data_dict[“gmsh:physical”][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={“tetra”: tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
cells=[(“triangle”, triangle_cells)],
cell_data={“name_to_read”:[triangle_data]})
meshio.write(“mesh.xdmf”, tetra_mesh)

meshio.write(“mf.xdmf”, triangle_mesh)

Scaled variables

E = 210000000000
nu = 0.3
mu, lambda_ = Constant(E/(2*(1 + nu))), Constant(Enu/((1 + nu)(1 - 2nu)))
#mu = 80769200000
rho = 7800
load = 0 #beachte Vorzeichen
F = load #Kraft auf Objekt (eventuell muss hier force per volume gewählt werden)
gamma = -9.81 #beachte Vorzeichen*
#beta = 121154000000
#lambda_ = beta
g = gamma

Create mesh and define function space

mesh = Mesh()
with XDMFFile(“mesh.xdmf”) as infile:
infile.read(mesh)
File(“winkel_mesh.pvd”).write(mesh)

mvc = MeshValueCollection(“size_t”, mesh, 1)
with XDMFFile(“mf.xdmf”) as infile:
infile.read(mvc, “name_to_read”)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File(“winkel_facets.pvd”).write(mf)

V = VectorFunctionSpace(mesh, ‘P’, 1)

Define boundary condition

bcl = DirichletBC(V, Constant((0, 0, 0)), mf, 1)
bcr = DirichletBC(V, Constant((0, 0, 0)), mf, 3)
bcup = DirichletBC(V, Constant((0, 200, 0)), mf, 2)

###bcup is displacement? how do I get a Force to this physical surface?

bcs = [bcl, bcr, bcup]

Define strain and stress

def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))

def sigma(u):
return lambda_nabla_div(u)Identity(d) + 2muepsilon(u)

Define variational problem

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, rho*g+F, 0))
T = Constant((0, 0, 0)) #traction force
#Z = Constant((0, 100, 0)) #force on physical plane
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds #+ dot(Z, v)*ds

Compute solution

u = Function(V)
solve(a == L, u, bcs)

Plot stress

s = sigma(u) - (1./3)*tr(sigma(u))Identity(d) # deviatoric stress
von_Mises = sqrt(3./2inner(s, s))
V = FunctionSpace(mesh, ‘P’, 1)
von_Mises = project(von_Mises, V)

Compute magnitude of displacement

u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print(‘min/max u:’,
u_magnitude.vector().get_local().min(),
u_magnitude.vector().get_local().max())

Save solution to file in VTK format

File(‘winkel/displacement_fenics.pvd’) << u
File(‘winkel/von_mises_fenics.pvd’) << von_Mises
File(‘winkel/magnitude_fenics.pvd’) << u_magnitude

[michi-ost]

Code for mesh: (Gmsh)

Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0.003, 0, 0, 1.0};
//+
Point(3) = {0.0165, 0, 0, 1.0};
//+
Point(4) = {0.0165, -0.003, 0, 1.0};
//+
Point(5) = {-0.0135, -0.003, 0, 1.0};
//+
Point(6) = {-0.0135, 0, 0, 1.0};
//+
Point(7) = {0, 0.015, 0, 1.0};
//+
Point(8) = {0.003, 0.015, 0, 1.0};
//+
Point(9) = {0.0275, 0.015, 0, 1.0};
//+
Point(10) = {0.0275, -0.015, 0, 1.0};
//+
Point(11) = {-0.0225, -0.015, 0, 1.0};
//+
Point(12) = {-0.0225, 0.015, 0, 1.0};
//+
Line(1) = {11, 10};
//+
Line(2) = {10, 9};
//+
Line(3) = {9, 8};
//+
Line(4) = {8, 2};
//+
Line(5) = {2, 3};
//+
Line(6) = {3, 4};
//+
Line(7) = {4, 5};
//+
Line(8) = {5, 6};
//+
Line(9) = {6, 1};
//+
Line(10) = {1, 7};
//+
Line(11) = {7, 12};
//+
Line(12) = {12, 11};
//+
Line(13) = {7, 8};
//+
Line Loop(1) = {12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
//+
Plane Surface(1) = {1};
//+
Line Loop(2) = {10, 13, 4, 5, 6, 7, 8, 9};
//+
Plane Surface(2) = {2};
//+
Extrude {0, 0, 0.020} {
Surface{2};
}
//+
Extrude {0, 0, 0.020} {
Surface{1};
}

//+
Extrude {0, 0.015, 0} {
Surface{4};
}
//+
Physical Surface(1) = {12};
//+
Physical Surface(2) = {29};
//+
Physical Surface(3) = {14};
//+
Physical Volume(4) = {1, 3};
//+
Physical Volume(5) = {2};

[michi-ost]

…I added the code to the onedrive link i’ve posted.

In addition I have an other question:

Where is fenics mostly used? Is it used in industry or is it more like a development program for research?
Is it uesd as a tool for simulation of complex geometry and much different parts and materials?

best regards

[dokken]

FEniCs is used for both research and industry. See for instance: https://arxiv.org/pdf/1804.10060.pdf
on large scale physical problems.
Im myself working on contact modell for FEniCS using multi-point constraints. I can let you know Once this is published.

For your case, it is easy to Connect the meshes,
Using boolean fragments in Gmsh. See:
Pygmsh Tutorial? for a simple use case.

If you Connect the meshes, the interface condition between the two objects will be that displacement on one side is equal to the displacement on the other side. (Both normal and tangential components)

Forces on surface is the same as applying traction on the top surface. This can be done using Measures and subdomains, see for instance:
Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio for a minimal code using a restricted surface measure.

[michi-ost]

Thanks for your advice.

1. Is this not similar to the code I wrote? (…similar to that what i did with the gmsh gui)

2. I applied the force with measures and subdomains like you said (…i tried, not really sure if I did it right)

But the result in paraview looks really strange…I have no idea why…

Define variational problem

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 200, 0)) #traction force
ds_custom = Measure(‘ds’, domain=mesh, subdomain_data=mf, subdomain_id=2)

a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

[dokken]

In your drawing, i did not see that there is space between the objects at one point. You also have to think about what interface conditions you would like inbetween the two domains. It could be:

1. They are full connected, i.e. u(domain_0)=u(domain_1) at the interfaces.
2. You allow slip, i.e. dot(u_domain_0, n)=(u_domain_1, n) at the interfaces (This would allow sliding in tangential directions). This condition is not simple to implement.

[Jesus_Amador]

@dokken

at there is space between the objects at one point. You also h

Hello Dokken,

I am also interested in learning how to implement in FEniCs a contact problem between two solids in which the contact zone is not known a priori.

I don’t know if you finally developed something about it that can be seen or if you could give me some advice to be able to implement the contact.

My goal is to develop a wear model and I have seen that there are phase field models and I don’t know up to what point they could be implemented to avoid updating the mesh or if it is easy in FEniCs to be updating the mesh in every step of time (I guess not).

[dokken]

A place to start is: https://comet-fenics.readthedocs.io/en/latest/demo/contact/penalty.html
And
Boundary conditions for an liquid sphere

[Forever]

Hi Dokken,

is your work published?
https://comet-fenics.readthedocs.io/en/latest/demo/contact/penalty.html the example here is too simple, for example there is no friction.

[dokken]

Some of my work (on MPC, i.e. No friction) can be found at: https://fenics2021.com/slides/dokken.pdf
and the source code is available at: GitHub - jorgensd/dolfinx_mpc: Extension for dolfinx to handle multi-point constraints.

I am currently working on a Nitsche approach, based on the work of Mlika et al.
but it is not ready for the public yet.

You could also have a look at:

evzenkorec/thesis_contact

work by @evzen on contact Mechanics with fenics.