How to correctly apply contact forces to a dynamic self-colliding mesh?

variational formulation
1

Hi,

I am implementing contact mechanics within a dynamic elastic model.
Two subparts of my mesh are colliding in the time.

So, in order to compute contact forces to avoid collision, I am applying the below process, but the contact forces do not seem to be taken into account.
Here are the main steps of the process:

  • I define two submeshes corresponding to the colliding subparts:
submesh1 = SubMesh(mesh, subdomains, 1)
submesh2 = SubMesh(mesh, subdomains, 2)
  • Then, I define VectorFunctionSpaces on these submeshes and compute contact forces on each submesh domain:
V1 = VectorFunctionSpace(submesh1, "CG", 1)
V2 = VectorFunctionSpace(submesh2, "CG", 1)

contactforces_submesh1 = Function(V1)
contactforces_submesh2 = Function(V2)

contactforces_submesh1.assign( compute_contact_force_1(time=0.) )
contactforces_submesh2.assign( compute_contact_force_2(time=0.) )
  • I then “convert” these forces from forces defined on V1 and V2 to forces defined on V (building a Vsub_2_V_dofmap), since the variational form I am using and associated unknown functions are defined on V:
V = VectorFunctionSpace(mesh, "CG", 1)
contactforces_submesh1_V = Function(V)
contactforces_submesh2_V = Function(V)

# convert contact forces in V1 space into contact forces in V space
for dofVsub in dofsV1:
   contactforces_submesh1_V.vector()[ V1_2_V_dofmap[dofVsub] ] = contactforces_submesh1.vector()[ dofVsub ] 

# convert contact forces in V2 space into contact forces in V space
for dofVsub in dofsV2:
   contactforces_submesh2_V.vector()[ V2_2_V_dofmap[dofVsub] ] = contactforces_submesh2.vector()[ dofVsub ]
  • Finally, I insert these contact forces, functions of V, into the definition of the variational form.
a(u,v) - L(v) - dot(contactforces_submesh1_V, v) * dx(1) + dot(contactforces_submesh2_V, v) * dx(2) = 0
# u = Function(V); v=TestFunction(V))
  • Then, at each time step, I update the contact forces functions of V1 and V2 only, but not the ones functions of V (since the second are supposed to depend on the first):
contactforces_submesh1.assign( compute_contact_force_1(time=t) )
contactforces_submesh2.assign( compute_contact_force_2(time=t) )
...
solver.solve()

Issue:
→ The problem is the two parts go on colliding in time… I manage to plot the contact forces both defined on V1, V2 and V, so I am supposing the updated forces are either not taken into consideration by the variational form, either not well updated in the time loop.

Could you please help me?

Many thanks,
Anne