Simple surface contact between different volumes

General
1

Hello all,

I have a very simple structure with three cylinders. This is the geo file:

//+
SetFactory("OpenCASCADE");
//+
Cylinder(1) = {0.0, 0, 0, 0, 0, 30, 0.5, 2*Pi};
//+
Cylinder(2) = {0.0, 1.0, 0, 0, 0, 30, 0.5, 2*Pi};
//+
Cylinder(3) = {0.8, 0, 0, 0, 0, 30, 0.3, 2*Pi};
//+
Physical Volume(10) = {2};
//+
Physical Volume(11) = {1};
//+
Physical Volume(12) = {3};
//+
Physical Surface(13) = {1};
//+
Physical Surface(14) = {2};
//+
Physical Surface(15) = {3};
//+
Physical Surface(16) = {4};
//+
Physical Surface(17) = {5};
//+
Physical Surface(18) = {6};
//+
Physical Surface(19) = {7};
//+
Physical Surface(20) = {8};
//+
Physical Surface(21) = {9};

I want to solve a simple elastic problem by fixing the top base of all cylinders (14, 17, 20) and apply a load along x direction in the bottom base of the smallest cylinder (21).
I use dolphin-converter to get the xml files.

This is the dolfin code:

from dolfin import *
from ufl import nabla_div
import numpy as np
import time

# Define the material properties of steel
E = 210e3   # Young's modulus
nu = 0.28   # Poisson's ratio

lame_mu = E/(2*(1+nu)) # Lame parameters
lame_lambda = (E*nu)/((1+nu)*(1-2*nu))

#############################################################

# Load mesh

mesh = Mesh('three_circle.geo.xml')
Volume = MeshFunction('size_t', mesh, "three_circle.geo_physical_region.xml")
bnd_mesh = MeshFunction('size_t', mesh, "three_circle.geo_facet_region.xml")

xdmf = XDMFFile(mesh.mpi_comm(),"Mesh.xdmf")
xdmf.write(mesh)
xdmf.write(Volume)
xdmf = XDMFFile(mesh.mpi_comm(),"boundaries.xdmf")
xdmf.write(bnd_mesh)
xdmf.close()
mesh= Mesh()
xdmf = XDMFFile(mesh.mpi_comm(), "Mesh.xdmf")
xdmf.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("Mesh.xdmf") as infile:
    infile.read(mvc, "f")

Volume = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc2 = MeshValueCollection("size_t", mesh, 2)

with XDMFFile("boundaries.xdmf") as infile:
    infile.read(mvc2, "f")

bnd_mesh = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = dx(metadata={'quadrature_degree': 2},
        subdomain_data=Volume)
ds = ds(metadata={'quadrature_degree': 2},
        subdomain_data=bnd_mesh)


##########################################################################



# Define function space
V = VectorFunctionSpace(mesh, 'Lagrange', 1) 


Area = 0.269959 


# 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 lame_lambda*nabla_div(u)*Identity(d) + 2*lame_mu*epsilon(u)


# Dirichelet BCs

bc1 = DirichletBC(V, Constant((0, 0, 0)), bnd_mesh, 14)
bc2 = DirichletBC(V, Constant((0, 0, 0)), bnd_mesh, 17)
bc3 = DirichletBC(V, Constant((0, 0, 0)), bnd_mesh, 20)

bcs =  [bc1, bc2, bc3]


# Boundary domain
ds = Measure('ds', subdomain_data=bnd_mesh)

load_force = 5

g_B= Constant((-load_force/Area, 0.0, 0.0))




# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(v, g_B)*ds(21)
start_time = time.time()

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

end_time = time.time()

elapsed_time = end_time - start_time

print("Elapsed time:", elapsed_time, "seconds")



# Compute von Mises stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = sqrt(3./2*inner(s, s))


V = FunctionSpace(mesh, 'Lagrange', 1)
von_Mises = project(von_Mises, V)



u.rename('Displacemente field vector (mm)', '')
von_Mises.rename('von Mises stress (MPa)', '')

cylinder_steel_file = XDMFFile('3_cylinders.xdmf')
cylinder_steel_file.parameters["flush_output"] = True
cylinder_steel_file.parameters["functions_share_mesh"] = True
cylinder_steel_file.write(u, 0.0)
cylinder_steel_file.write(von_Mises, 0.0)

This code is working but I can see that the smallest cylinder is penetrating one of the other cylinders. What would be the best approach to avoid this? I have done my research and found that I could include a contact penalty. Is this feasible for problems with many physical volumes (ex: 10 cylinders touching each other).

2

Contact is not a simple problem if you want it to work with MPI.
Some references can be found at: How can I correct collisions and preserve mesh during dynamic simulation? - #8 by dokken

3

Thank you for your answer. I see that there are more demos for dolfinx about contact mechanics. Would you suggest that I switch from dolfin to dolfinx?

4

It depends on your familiarity with DOLFIN.
I would look at: The new DOLFINx solver is now recommended over DOLFIN
and look at the actual resources to see what will fit you best.

DOLFIN is no longer maintained (in the sense that there are no-one fixing bugs, and rarely anyone making PRs to ensure that it is compatible with new compilers etc).

5

Would it be feasible to implement a penalty between different surfaces?
In this case I know that surface 13 will collide with try to penetrate surface 19.
If I have many surfaces in which I know beforehand that they will have some kind of contact, would it be reasonable to include a penalty for all of them?

6

As long as you are able to compute the distance between those surfaces, and put those into a dolfin.Function it should be feasible.

7

Thanks! I will give a look at all demos with contacts and see what I can do.