Yade

Unexpected PFacet Collision Physics

Asked by Jeffrey Knowles

I am using Yade 2020.01a. The MWE below simulates the collision of two PFacets. After PFacet 1 (on the left) collides with PFacet 2 (on the right) unrealistic physical behavior is observed. I would expect PFacet 2 to rotate clockwise about the y-axis. However, it appears that PFacet 2 has a tendency to rotate counter-clockwise about the z-axis with some angular momentum about the y-axis (but in the counter-clockwise direction). I have tried smaller and smaller time steps, but observe the same behavior.

I am not sure if this is a bug, or if I have a misconception about PFacets. At any rate, I would appreciate any time that could be allocated in addressing this issue.

MWE:

import sys
p = os.getcwd() # current working directory
sys.path.append(p)
from yade import pack, geom, PyRunner, qt
from yade.gridpfacet import *

############################################# engines ###########################################
O.engines=[

 ForceResetter(),
 InsertionSortCollider(

  [
  Bo1_PFacet_Aabb(),
  Bo1_GridConnection_Aabb()
  ]),

 InteractionLoop(

  [
  Ig2_GridConnection_PFacet_ScGeom(),
  Ig2_Sphere_PFacet_ScGridCoGeom(),
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
  Ig2_PFacet_PFacet_ScGeom()
  ],

  [
  Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  ],

  [
  Law2_ScGeom_FrictPhys_CundallStrack(),
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm = True),
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
  Law2_ScGridCoGeom_CohFrictPhys_CundallStrack(),
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
  ]),

  NewtonIntegrator(gravity=(0,0,0), damping=0.0)
  ]

O.dt = 1.0e-9 # times step

############################################# materials ###########################################
cmatName = 'cmat'
fmatName = 'fmat'

O.materials.append( CohFrictMat( young=1.0e14,poisson=0.3,density=1000.0, normalCohesion=1.0e12 ,shearCohesion=1.0e12, fragile = False, momentRotationLaw=True, isCohesive = True, etaRoll=1000.0, etaTwist = 1000.0, frictionAngle=0.0, alphaKr = 100.0, alphaKtw = 100.0, label=cmatName ) )

O.materials.append( FrictMat( young=1.0e14,poisson=0.3,density=1000.0,frictionAngle=0.0,label=fmatName ) )

############################################## pfacets #############################################

# pfacet 1
gnh1_0 = O.bodies.append( gridNode([0.0, 0.0, 0.3],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh1_2 = O.bodies.append( gridNode([-1.0, 1.0, 0.3],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh1_3 = O.bodies.append( gridNode([-1.0, -1.0, 0.3],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )

O.bodies.append( gridConnection(gnh1_2,gnh1_3,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh1_0,gnh1_2,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh1_0,gnh1_3,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )

O.bodies.append( pfacet(gnh1_0, gnh1_2, gnh1_3,material=fmatName,color=[1, 0, 0]) )

# assign initial velocity to pfacet 1
vi = Vector3(10.0, 0.0, 0.0)
for b in O.bodies:
 b.state.vel = vi

# pfacet 2
gnh2_0 = O.bodies.append( gridNode([0.21, 0.0, 0.0],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh2_1 = O.bodies.append( gridNode([0.21, 1.0, 1.0],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh2_2 = O.bodies.append( gridNode([0.21, -1.0, 1.0],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )

O.bodies.append( gridConnection(gnh2_0,gnh2_1,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh2_0,gnh2_2,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh2_1,gnh2_2,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )

O.bodies.append( pfacet(gnh2_0, gnh2_1, gnh2_2,material=fmatName,color=[1, 0, 0]) )

#################################### appearance and view ##########################################
qtr = qt.Renderer()
qt.Controller()
v = qt.View()

qtr.bgColor = [1,1,1]

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Launchpad Janitor (janitor) said : 		#1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Revision history for this message
Jeffrey Knowles (knowleje) said : 		#2

Hello,

I would like to reopen this question.

Thank you,

Jeff Knowles

On Fri, Jan 1, 2021 at 5:06 AM Launchpad Janitor <
<email address hidden>> wrote:

> [This email originated from outside of OSU. Use caution with links and
> attachments.]
>
> Your question #694548 on Yade changed:
> https://answers.launchpad.net/yade/+question/694548
>
> Status: Open => Expired
>
> Launchpad Janitor expired the question:
> This question was expired because it remained in the 'Open' state
> without activity for the last 15 days.
>
> --
> If you're still having this problem, you can reopen your question either
> by replying to this email or by going to the following page and
> entering more information about your problem:
> https://answers.launchpad.net/yade/+question/694548
>
> You received this question notification because you asked the question.
>

Revision history for this message
Janek Kozicki (cosurgi) said : 		#3

This seems related to:

https://answers.launchpad.net/yade/+question/692998
https://gitlab.com/yade-dev/trunk/-/issues/167

but not sure if that's really related to it.

Revision history for this message
Jeffrey Knowles (knowleje) said : 		#4

Hello Janek,

Thank you for your reply. I do not believe the links you sent are really
that related to the issue.

When I run a simulation where a grid node of the left body collides with a
grid node of the right body, the expected angular velocities are observed.
However, if the grid node of the left body collides with the face of the
right body, then the expected angular velocities are not observed (not even
qualitatively).

I think it may be due to the barycentric distribution of force in equations
(15) and (17) of [1]. Please let me know your thoughts.

[1]
https://www.sciencedirect.com/science/article/pii/S0266114415001235?casa_token=maargb3eHr4AAAAA:LH_-_RVKMaiC7srcr5hL2J_cziD14z0mWc3ovLpPlYhZv8n27RaFA9E_fh-_Uj-41VBCC_r5s-4

Thank you,

Jeff

On Mon, Jan 4, 2021 at 7:16 AM Janek Kozicki <
<email address hidden>> wrote:

> [This email originated from outside of OSU. Use caution with links and
> attachments.]
>
> Your question #694548 on Yade changed:
> https://answers.launchpad.net/yade/+question/694548
>
> Janek Kozicki posted a new comment:
> This seems related to:
>
> https://answers.launchpad.net/yade/+question/692998
> https://gitlab.com/yade-dev/trunk/-/issues/167
>
> but not sure if that's really related to it.
>
> --
> You received this question notification because you asked the question.
>

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#5

Hi,

Sorry for not checking your question earlier. I'll try and check the script.

Could you elaborate the comment on barycentric distribution? Did you
spot an error in some equation?

Bruno

Revision history for this message
Jeffrey Knowles (knowleje) said : 		#6

Hello Bruno,

Not a problem.

I am not sure if there is an error in the barycentric equation. Rather, I
think it is using the barycentric equation itself which creates unphysical
results. The angular momentum is not conserved during a node-on-face
collision because, if I understand correctly, the force is not applied to
the collision point on the face, but is instead distributed to the nodes of
the pfacet. Therefore, when summing torques, the moment arms are crossed
with the forces at the nodes, instead of the collision point on the face.
So, the total force is the same (equations 15 and 17 of [1]) when summed
over the nodes, but the torques are different than when calculating them in
the usual sense (r_face x F not necessarily equal to sum(r_node_i x F_i) ).
Therefore, it might not be a bug or typo, but maybe I need to modify the
source code so forces are not distributed to the nodes. If so, could you
direct me to the relevant files?

Thank you,

Jeff

On Tue, Jan 5, 2021 at 6:35 AM Chareyre <
<email address hidden>> wrote:

> [This email originated from outside of OSU. Use caution with links and
> attachments.]
>
> Your question #694548 on Yade changed:
> https://answers.launchpad.net/yade/+question/694548
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
> Hi,
>
> Sorry for not checking your question earlier. I'll try and check the
> script.
>
> Could you elaborate the comment on barycentric distribution? Did you
> spot an error in some equation?
>
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/694548/+confirm?answer_id=4
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/694548
>
> You received this question notification because you asked the question.
>

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#7

I would say the angular momentum is (should be...) conserved.
Simply, the PFacet is, from an inertial point of view, a triangle with a solid sphere at each vertex. As if all the rest of the volume was occupied by a solid of mass density ~0. This is peculiar but not unphysical in itself, it shouldn't give something strange by design.

Defining a unique torque using the single contact point (i.e. not distributing) is possible only if the facet is a rigid clump, not your case. You have a deformable facet with 3x3 rotational DOFs. Then you need 3 torques.

> r_face x F not necessarily equal to sum(r_node_i x F_i)

I see what you mean... but they _should_ be necessarily equal if we put the equations correctly. I'll have to go back to that paper. I suggest you also check the source, since it coud differ from paper.
Most of the code is in common/PFacet.hpp,cpp

HTH
Bruno

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#8

Tries your script. To me it looks like there is a wrong numbering of nodes, a mathematical error in the force/torque distribution, or something alike. Distributing is not the problem in itself. I'll have a look.
B

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#9

At impact time the forces are the same on nodes 7 and 8 while based on the symmetry it should be the same on 8 and 9 instead.
Really looks like at some point the ordering of the nodes is failed. Now it needs to find how.

In [1]: O.forces.f(7)
Out[1]: Vector3(266280.2931349052233,23.39075633919068764,-20.63102091593614418)

In [2]: O.forces.f(8)
Out[2]: Vector3(266198.2480056955828,23.58386741905894368,-21.62634969041041444)

In [3]: O.forces.f(9)
Out[3]: Vector3(1242204.077827499015,115.5724137135645435,-101.9188704456528143)

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#10

issue on gitlab for tracking progress on the bug: https://gitlab.com/yade-dev/trunk/-/issues/184

Let me know if you make progress in understanding the issue. I think I may fix it but I have no ETA, I hope it is not a blocker for you.

Bruno

Revision history for this message
Jeffrey Knowles (knowleje) said : 		#11

Hello Bruno,

Thank you for all of your time and expertise on this issue.

I have run some full-scale models with these pfacet bodies (essentially a
large packing of these bodies in a wire cage), and found that the linear
momentum is significantly more important than the angular momentum when
calibrating my model with laboratory results, so it is not the end of the
world.

At any rate, I would be interested to hear if you fix it.

Best,

Jeff

On Mon, Jan 11, 2021 at 11:51 PM Bruno Chareyre <
<email address hidden>> wrote:

> [This email originated from outside of OSU. Use caution with links and
> attachments.]
>
> Your question #694548 on Yade changed:
> https://answers.launchpad.net/yade/+question/694548
>
> Bruno Chareyre posted a new comment:
> issue on gitlab for tracking progress on the bug: https://gitlab.com
> /yade-dev/trunk/-/issues/184
>
> Let me know if you make progress in understanding the issue. I think I
> may fix it but I have no ETA, I hope it is not a blocker for you.
>
> Bruno
>
> --
> You received this question notification because you asked the question.
>

Revision history for this message
Launchpad Janitor (janitor) said : 		#12

This question was expired because it remained in the 'Open' state without activity for the last 15 days.