sphere-facet contact

Asked by ytang

Hi all,

As we all know, the timestep is calculated by p-wave analysis for the sphere-sphere/sphere-facet contact[1].

The last time, I posted a qeustion[2], I said that the reason for the particle blows up is related to the timestep. when the number of the penetrator segments is 10, the simulation runs well, while when the number of the penetrator segments is 50, the particles blow up. (the other parameters are the same, except for the number of the segments). This means the timestep for these two simulations should be the same. why the particles blow up?
######################################################
I still didn't fully understand the situation of the sphere-facet contact.
######################################################
Here, I just want to express my understanding and confusion about this question.
(1) If a sphere contact with the facet, there is one contact point if the central point of the sphere is located within the triangular area of the facet.
(2) There are two contact points if the sphere contact at the edges of two neighboring facets (In convex edge). At the same time, there are two contact forces, one is from the sphere-facet1, the other is from sphere-facet2.
(3) Sometimes we said that the contact normal is perpendicular to the facet. the precondition for this is that the contact is within the triangular area of the facet. If the contact point is between the sphere and the edge (which is made by two facets), the contact normal will not be perpendicular to the facet.
(4) As for the vertices, what are the differences if the cone is made of 10 segments and 50 segments? (the higher the segments, the higher the resultant force on the vertice???)
(5) Jan mentioned that it is like the cone tip has 50x higher stiffness. But actually, the stiffness is still the same, right??? If this is correct. how do we determine the timestep for the sphere-facet contact?

overall, there might be two or more contacts between the sphere-edge/sphere-vertice contact.
###################################################################
I also read some articles about pfacet[3]. Here, it mentioned that there is only one contact between the sphere and the edge (convex edge).
############################################
Do we have any detailed information about the sphere-facet contact? where I can understand it thoroughly.

References:
[1]https://yade-dem.org/doc/yade.utils.html?highlight=pwave#yade._utils.PWaveTimeStep
[2]https://answers.launchpad.net/yade/+question/694208
[3]https://www.sciencedirect.com/science/article/abs/pii/S0266114415001235

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said :
#1

Hello,

> This means the timestep for these two simulations should be the same. why the particles blow up?

It depends what "should be the same" does mean.
Yes, the time step is ok for sphere-one facet contact and should be the same.
You have sphere-multiple facet contact, so the overall time step should NOT be he same, see below.

All the points you mentioned are basically correct.

The stiffness of sphere-one facet contact has absolutely no problems.
Contact point, normal vector and forces are computed "naturally" and expectedly.
Also critical time step for the contact follows the p-wave rules.

What is a problem is contact of one sphere with multiple facets at once, where in reality there is just one surface.
Time step computation is derived from mass-spring system, from **one** interaction.

At the edge, "normally" the sphere would interact with the surface with its stiffness K.
For penetration depth pd, the repulsive force F would be something like F=K*pd.
But you have two facets. The sphere get repulsive force from both facets independently, from each one F1=K*pd, totally from two facets is F=2*F1=2*K*pd.

Now if you have vertex with 50 facets, the sphere get repulsive force from 50 facets independently, from each one F1=K*pd, totally from 50 facets is F=50*F1=50*K*pd.

F1 = K*pd
F = 50*K*pd = (50*K)*pd
K, true facet stiffness is still the same.
But the "overall vertex stiffness" is (50*K), 50 times higher, it should be clear from the formulas why.

So you have discrepancy between critical time step assumptions (one interaction) and reality (50 interactions).

> Do we have any detailed information about the sphere-facet contact? where I can understand it thoroughly.

Source code is probably the best source. But as mentioned, the sphere-facet contact itself is not the problem.
The problem is actually not related to sphere-facet contact, but to the problem of contact of one sphere with multiple bodies "at the same place".
If you model some rough surface with multiple extensively overlapping spheres, you would come to the same problem.

Also, in 3D, spheres may have multiple interactions in "random" directions, resulting on "overall stiffness" higher than that of one interaction. Therefore safety factor is often used (with value of 0.7 or 0.8) w.r.t PWaveTimeStep result.

cheers
Jan

Revision history for this message
Jérôme Duriez (jduriez) said :
#2

Hi,

> As we all know, the timestep is calculated by p-wave analysis for the sphere-sphere/sphere-facet contact

I actually do not know this but rather rely (for spherical particles at least) on GlobalStiffnessTimeStepper for the "best" (greatest, and still less than critical) estimate of time step.

The corresponding algorithm directly involves properties of the interaction network, which does control critical time step alongside particles' masses and FrictMat.young, as you evidenced here with your case.

See eg https://yade-dem.org/doc/formulation.html#general-mass-spring-system and the following paragraph

Revision history for this message
Jan Stránský (honzik) said :
#3

@Jerome: I didn't know that GlobalStiffnessTimeStepper considers also directions and accounts for this case, thanks for edification!

@Ytang: Even if GlobalStiffnessTimeStepper solves the time step problem, it probably slows down significantly the whole simulation just because how the tip is modelled.. I think PFacets are worth trying.

cheers
Jan

Revision history for this message
ytang (ytang116) said :
#4

Thanks Jan Stránský, that solved my question.

Revision history for this message
ytang (ytang116) said :
#5

Hi Jan and Jerome,

thank you!
ytang