# Detecting collisions in FEM

Say I want to study the deformation of a pitchfork when you have it fixed on the bottom and push one side.

<< NDSolveFEM
Ω =
RegionDifference[Cuboid[{0, -5}, {1, 0}],
Cuboid[{0.45, -4.5}, {.55, 0}]];
Ω // DiscretizeRegion
bcs = DirichletCondition[#, {z <= -4.99}] & /@ {u[y, z] == 0.,
v[y, z] == 0.};
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 0.005]
planeStress = {Inactive[
Div][{{0, -((Y*ν)/(1 - ν^2))}, {-(Y*(1 - ν))/(2*(1 \
- ν^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y/(1 - ν^2)),
0}, {0, -(Y*(1 - ν))/(2*(1 - ν^2))}}.Inactive[Grad][
u[y, z], {y, z}], {y, z}],
Inactive[
Div][{{0, -(Y*(1 - ν))/(2*(1 - ν^2))}, {-((Y*ν)/(1 \
- ν^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y*(1 - ν))/(2*(1 - ν^2)),
v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, ν -> 33/100};


That's my domain

and this my results:

{uif, vif} =
NDSolveValue[{planeStress == {NeumannValue[1, y == 0 && z > -.1],
0}, DirichletCondition[u[y, z] == 0, z == -5],
DirichletCondition[v[y, z] == 0, z == -5]}, {u,
v}, {y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


I was expecting that the right arm of the pitchfork would somehow know when the other one was touching it, but this is not the case at all.

I don't have much experience with FEM so I don't know if what I'm asking is impossible due to the nonlocality of the problem or what not. Is there a way to obtain the correct behaviour in mma, ie how to make the left arm collide with the right one and that they both deform due to the force applied only on the left one?

• The topology of the solution region changes when the two arms come into contact. You may have to terminate the calculation (use WhenEvent when the this event occurs, then use the output as initial condition for a continued calculation with the two arms joined. Even this may not be quite correct, because the contact slides. Interesting question. – bbgodfrey Aug 19 '16 at 18:27
• mmmh, one could integrate till it touches for the first time, then compute the force it applies on the right arm, move the right arm, then integrate the left arm till the next contact and proceed like that till the balance of forces is reached? That could be an option to avoid resolving the time dependent problem... – tsuresuregusa Aug 19 '16 at 20:19
• The equations you have given to the Mathematica solver to discretize and solve via F.E.M don't take into account collisions/overlapping in the mesh. Adding collisions must be done via crafting boundary conditions or penalty forces which change dynamically. This is definitely non-trivial: it's quite complicated, and is normally done iteratively. Additional note: Furthermore the Mathematica solver linearizes your equations. The linearized equation is only valid in the small deformation case: distances and angles barely change (imagine for instance 10 m steel beam bending a few millimeters) ... – darmual Aug 7 '18 at 11:50
• [second part of the comment]: ... which is definitely not your case. So your results, even without the overlapping, wouldn't be valid. You could see a clear example of this by changing the Neumann value to -50. – Kuba Aug 7 '18 at 13:09
• @darmual I moved your answer to comments section since it didn't really provide an answer – Kuba Aug 7 '18 at 13:10

Here it is necessary to determine not the impact, but the additional stress arising at the point of contact of two elastic elements. The simplest solution method is to determine the force that acts on each element in contact. I will show the simplest code that can be complicated to infinity.

<< NDSolveFEM
Ω = RegionDifference[Cuboid[{0, -5}, {1, 0}], Cuboid[{0.45, -4.5}, {.55, 0}]];

Ω // DiscretizeRegion

bcs = DirichletCondition[#, {z <= -4.99}] & /@ {u[y, z] == 0., v[y, z] == 0.};
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 0.005];
planeStress = {Inactive[
Div][{{0, -((Y*ν)/(1 - ν^2))}, {-(Y*(1 - ν))/(2*(1 \
- ν^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y/(1 - ν^2)),
0}, {0, -(Y*(1 - ν))/(2*(1 - ν^2))}}.Inactive[Grad][
u[y, z], {y, z}], {y, z}],
Inactive[
Div][{{0, -(Y*(1 - ν))/(2*(1 - ν^2))}, {-((Y*ν)/(1 \
- ν^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y*(1 - ν))/(2*(1 - ν^2)),
v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, ν -> 33/100};
{uif, vif} =
NDSolveValue[{planeStress == {NeumannValue[1, y == 0 && z > -.1] +
NeumannValue[-1/3, y == 0.45 && z > -.1] +
NeumannValue[1/3, y == 0.55 && z > -.1], 0},
DirichletCondition[u[y, z] == 0, z == -5],
DirichletCondition[v[y, z] == 0, z == -5]}, {u,
v}, {y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}


Note that in this solution, the elements do not crawl together, but are deformed together. Consider the general case of an arbitrary combination of parameters of force and elasticity.

<< NDSolveFEM
\[CapitalOmega] =
RegionDifference[Cuboid[{0, -5}, {1, 0}],
Cuboid[{0.45, -4.5}, {.55, 0}]];
\[CapitalOmega] // DiscretizeRegion;
bcs = DirichletCondition[#, {z <= -4.99}] & /@ {u[y, z] == 0.,
v[y, z] == 0.};
mesh = ToElementMesh[\[CapitalOmega],
"MaxCellMeasure" -> 0.005]; f = 2;
planeStress = {Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
u[y, z], {y, z}], {y, z}],
Inactive[
Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, {-((Y*\[Nu])/(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)),
v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, \[Nu] -> 33/100};
sol = ParametricNDSolveValue[{planeStress == {NeumannValue[f,
y == 0 && z > -.1] + NeumannValue[-g, y == 0.45 && z > -.1] +
NeumannValue[g, y == 0.55 && z > -.1], 0},
DirichletCondition[u[y, z] == 0, z == -5],
DirichletCondition[v[y, z] == 0, z == -5]}, {u,
v}, {y, z} \[Element] mesh, {g}];

sol1 =
ParametricNDSolveValue[{planeStress == {NeumannValue[f,
y == 0 && z > -.1] + NeumannValue[-g, y == 0.45 && z > -.1] +
NeumannValue[g, y == 0.55 && z > -.1], 0},
DirichletCondition[u[y, z] == 0, z == -5],
DirichletCondition[v[y, z] == 0, z == -5]},
u, {y, z} \[Element] mesh, {g}];

g0=g/.FindRoot[
sol1[g][.45, 0] - sol1[g][.55, 0] == .1, {g, f/2}] // Quiet

Out[]= 0.834936
dmesh = ElementMeshDeformation[mesh, sol[g0],
"ScalingFactor" -> 1];
{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}


For f=1, we find g0=0.335075, which is close to 1/3 found by another method.

• Could you explain how you did "determine the force that acts on each element in contact"? – user21 Nov 20 '18 at 9:39
• @user21 I took it as a first approximation. You can further complicate the code using contact theory. – Alex Trounev Nov 20 '18 at 12:55
• How did you take this first approximation? – user21 Nov 20 '18 at 13:20
• @user21 I searched for the minimum value when there is no creep and there is contact. 1/3 is close to the exact solution in this problem. – Alex Trounev Nov 20 '18 at 13:41
• This is much better. Thank you. Now it's clear how your method works. – user21 Nov 21 '18 at 7:02