# PFacet model - contact data mining

Hello there!

Recently I've been working a lot with the PFacet approach to obtain an elastic model for bigger but simple bodies in DEM simulations.

I managed to obtain correct mechanical behaviour of my elastic body. Now I am researching the contact with external particles.

I made up a simple example of a PFacet cuboid that is penetrated by a sphere. The goal is to do that at every different position on the cuboid (edge, corner, intersection...) to verify that at all positions the same contact behaviour (penetration depth, force etc.) appears. I do this in order to get a better understanding of how exactly the model works (coding wise) as I want to change the contact model to an adapted hertz-mindlin model afterwards.

When I perform this penetration test e.g. within the area of a PFacet and on the intersection of 2 PFacets (so directly on a cylinder connection) I can see that the behaviour is according to the Cundall Strack law, which was expected, and that I only have 1 interaction on the sphere.

However, if I perform this test on the very edge of the cuboid, like the following script shows, I get unexpected results.

Mostly there is no interaction with the sphere, but according to the position and velocity of the sphere it behaves correctly and identical to other sphere positions (as mentioned before).

If I go through the simulation step by step, at some timesteps 3 interactions are found. Those 3 interactions are between the sphere and 3 different PFacets, one of them should be not even close to the sphere. If 3 interactions are found, I would have expected that there are 2 interactions with PFacets and 1 with the cylindrical connection.

Still, this isn't in accordance with the other tests where the sphere is at a different position and I am still not able to extract the wanted data.

So my question is: How can I correctly mine the contact parameters like penetration depth, force on the sphere and number of interactions for my penetration tests?

Am I using the wrong functions or am I looking at the wrong spot?

Thanks in advance!

P.S.: Since it will be needed in the future: If there is an "easy" way (only changing python code, not c++ source code) to change the contact law for the PFacet - external particles interactions, I would be happy to receive hints.

Working with YADE 2020.01a on Ubuntu 20.04.2 LTS

#######

from yade import plot, qt

import os, sys, time

from yade.gridpfacet import *

import numpy as np, math

""" All manual entries for the simulation """

O.dt = 1e-8 # timestep

node_mat_name = 'gridNodeMat'

node_r = 10 * 1e-3 # radius for the nodes, cylinders and pfacets

node_E = 52.*1e9 # (concrete)

node_poisson = 0.167

node_rho = 2750

node_phi = 35.

node_normal_coh = 3e1000 # high values to not lose node interactions. total elasticity, no plasticity

node_shear_coh = 3e1000

pfacet_mat_name = 'pFacetMat' # this material has no influence on the beam stiffness, only for contacts

pfacet_E = node_E

pfacet_poisson = node_poisson

pfacet_rho = node_rho

pfacet_phi = node_phi

color = [0, 0, 1] # color for certain bodies (PFacets)

sphere_mat_name = 'sphereMat'

sphere_poisson = 0.2

sphere_phi = 24.

sphere_E = 30.*1e6

sphere_rho = 2660

sphere_r = 40 * 1e-3 # radius of the penetrationg sphere

sphere_v = -0.1*100 # velocity of the penetrating sphere

dim_x = 0.5 # length in x-direction of the cuboid

dim_y = 0.6 # width of the cuboid

dim_z = 0.7 # height of the cuboid

# initial position of the penetrating sphere

sphere_x = (1/2) * dim_x

sphere_y = 0

sphere_z = dim_z + sphere_r + node_r

engine_gravity = (0,0,0) # gravity acceleration

engine_damping = 0.0 # general damping factor

plot_label = "plotter"

plotData_period = 50 # iterations

"""" The simulation engine """

O.engines=[

ForceResett

InsertionSo

Interaction

),

NewtonInteg

PyRunner( command = 'addPlotData()', iterPeriod = plotData_period),

]

""" Materials """

O.materials.append(

CohFrictMat( young = node_E, poisson = node_poisson, density = node_rho,

label = node_mat_name

)

)

O.materials.append(

FrictMat(

young = pfacet_E, poisson = pfacet_poisson, density = pfacet_rho,

label = pfacet_mat_name

)

)

O.materials.append(

FrictMat(

young = sphere_E , poisson = sphere_poisson, density = sphere_rho,

label = sphere_mat_name

)

)

""" Bodies """

# adding the 8 corner nodes of the cuboid

nodesIds = [] # holds all gridnode ids

nodesIds.append( O.bodies.append( gridNode( [0,0,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,0,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [0,dim_y,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,dim_y,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [0,0,dim_z], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,0,dim_z], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [0,dim_y,dim_z], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,

# adding the connections and pfacets with the pfacetCreator3 function

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

# the penetrating sphere

idSphere = O.bodies.append( sphere( [sphere_

""" Movements """

O.bodies[

""" PyRunner Functions """

def addPlotData() :

# if len( O.bodies[

# penetrationDepth = 0

# force_ball = 0

# if O.iter > 200000 :

# O.save(workingDir + '/' + sim_save_folder + filename + '.yade.bz2')

# plot.saveGnuplo

# O.pause()

# else :

# penetrationDepth = O.bodies[

# force_ball = O.bodies[

plot.addData(

time = O.time,

nrInt = len( O.bodies[

# penetrationDepth = penetrationDepth,

# force_ball = force_ball,

pos_ball = O.bodies[

vel_ball = O.bodies[

forceZ_ball = O.forces.

)

""" Plots """

# plot.plots = {'penetrationDepth' : 'force_ball', 'time' : 'penetrationDepth', 'time ' : 'force_ball', 'time ' : 'nrInt'}

plot.plots = {'time' : 'nrInt', 'time ' : 'pos_ball', 'time ' : 'vel_ball', 'time ' : 'forceZ_ball'}

plot.plot( subPlots = True )

""" Viewer """

qt.Controller()

qtv = qt.View()

qtv.ortho = True

qtv.eyePosition = [0., 1.5, 0]

qtv.viewDir = [0, -1, 0]

qtv.axes = False

qtv.grid = (False, False, False)

qtr = qt.Renderer()

qtr.light2=True

qtr.lightPos=

qtr.bgColor=[1,1,1]

## Question information

- Language:
- English Edit question

- Status:
- Answered

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Last query:

- Last reply:

Jérôme Duriez (jduriez) said : | #3 |

Hi

> How can I correctly mine the contact parameters like penetration depth, force on the sphere and number of interactions for my penetration tests?

for cont being some contact interaction (eg cont = O.interactions[

See also https:/

> Since it will be needed in the future: If there is an "easy" way (only changing python code, not c++ source code) to change the contact law

I do not see one that would properly fit YADE design. DEM operations in YADE are executed from the C++ side, through C++ source of Engines and LawFunctors, for instance. Python is "just" for designing and interrogating what's going on, and that is already something marvelous.

You could still always rewrite all DEM operations in pure Python but that would be slower.

Paul Pircher (chabs300) said : | #4 |

Hi and thanks for your input!

I know how I am supposed to receive parameters of interactions.

In my case, as can be seen in the code in the "PyRunner Function" - section, I am using len( O.bodies[

Furthermore I am extracting the phys-norm-force from the first interaction with the sphere via O.bodies[

I also had a look at O.forces.

However, the problem that I am facing in this particular example with PFacets is, that using this approach does NOT provide correct data that are in accordance with the visual observations. NO INTERACTION can be found in here. Which is why I do not receive a penetrationDepth parameter.

Since I am using the same approach and functions as you mentioned, I assume that I the data mining itself would be correct, but my data is found somewhere else. And I don't know where.

Since the observation of position and velocity does show correct behaviour, somewhere the contact has to be computed.

Although there seems to be no interaction for that.

I hope that this is understandable.

> I do not see one that would properly fit YADE design.

Alright, thank you. This is also what I concluded. Seems like I am going to have a look at the pfacet source code in order to know what I need to adapt or add to obtain the behaviour of a different contact model.

Since the pfacet approach works fine for me so far (except that data mining issue mentioned above) I think, that only few adaptations need to be done to change the physical computations of the contact itself.

Thanks!

Paul Pircher (chabs300) said : | #5 |

Since I just had a look at more positions:

It seems like that only the contact of the sphere at the position of a gridconnection at the edge of the cuboid is an issue.

Depending from what direction the sphere penetrates the cuboid I either get zero or three interactions.

At any other position on the cuboid (somewhere in the middle of a side of the cuboid or even at the corner) I receive the expected 1 interaction.

I am still clueless about that edge contact though.

With zero interactions I don't know where my wanted data is stored (e.g. penetration depth).

With 3 interactions I don't know what interaction is the "right" one or how those multiple but similar contacts are handled.

As always: The contact behaviour according to position and velocity of the sphere works perfectly fine at every location.

Cheers!

Launchpad Janitor (janitor) said : | #6 |

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

Paul Pircher (chabs300) said : | #7 |

Can someone give me more insight here? :)

Hi Paul,

I did not check all details but here are quick answers. I would say the excess interactions you see are just virtual interactions. Try filtering them out by checking "interaction.

The zero-interaction case sounds like a real issue though. Probably linked to a couple issues we identified recently. I hope we can fix them soon. I may use your example to try and understand the problem.

> change the PFacet - external particles interactions

From the top of my head: if you assign another material to the PFacets and to the spheres, and if you use the corresponding functors (just like for changing the sphere-sphere contact model) then it should change the contact model.

Regards

Bruno

Paul Pircher (chabs300) said : | #9 |

Hi Bruno,

thanks for your answer!

Only a short response from my side as well.

The last time I checked all 3 interactions were real at the contact.

The corresponding bodies in contact with the penetrating sphere were only pfacets. One of the pfacet wasn't even close to the contact location. No interaction with a grid-connection that might be expected.

I found out that the case for zero interactions is more complex.

At even iteration numbers there were 0 interactions, at odd iteration numbers there were 3 interactions.

I noticed the post where general PFacet problems are discussed and had a look at it.

The thing is that I am often not sure if I found something of serious interest (bug, unintended feature, etc.) or my limited knowledge is the issue. I am not a software engineer and therefore I struggle with too complex code especially with C++.

So I rather post general questions here.

I will prepare my examples, get some screenshots and post them here again. Maybe that helps.

> From the top of my head: if you assign another material to the PFacets and to the spheres, and if you use the corresponding functors (just like for changing the sphere-sphere contact model) then it should change the contact model.

I need to try that again. If this works it would save me a ton of work.

Although I am using a contact model that was developed at my departement. So it isn't implemented in YADE straight away and therefore might cause some issues.

Best regards,

Paul

Paul Pircher (chabs300) said : | #10 |

Hi there,

now I got some details and screenshots ( https:/

I was right: All interactions were real. I plotted O.interactions.

18 for no external interaction is clear as the cube I used should have 18 interactions since each gridConnection enables one gridNode-gridNode interaction.

Hence, the sphere contact leads to 3 interactions at that edge position.

Also as seen on the screenshots, the interacting bodies are only PFacets with one not even close to the contact.

I only did very few changes in the script but for completeness in paste the script I used to generate these screenshots. Please note the iteration period for the plot is an odd number.

#######

from yade import plot, qt

import os, sys, time

from yade.gridpfacet import *

import numpy as np, math

""" All manual entries for the simulation """

O.dt = 1e-8 # timestep

node_mat_name = 'gridNodeMat'

node_r = 10 * 1e-3 # radius for the nodes, cylinders and pfacets

node_E = 52.*1e9 # (concrete)

node_poisson = 0.167

node_rho = 2750

node_phi = 35.

node_normal_coh = 3e1000 # high values to not lose node interactions. total elasticity, no plasticity

node_shear_coh = 3e1000

pfacet_mat_name = 'pFacetMat' # this material has no influence on the beam stiffness, only for contacts

pfacet_E = node_E

pfacet_poisson = node_poisson

pfacet_rho = node_rho

pfacet_phi = node_phi

color = [0, 0, 1] # color for certain bodies (PFacets)

sphere_mat_name = 'sphereMat'

sphere_poisson = 0.2

sphere_phi = 24.

sphere_E = 30.*1e6

sphere_rho = 2660

sphere_r = 40 * 1e-3 # radius of the penetrationg sphere

sphere_v = -0.1*100 # velocity of the penetrating sphere

dim_x = 0.5 # length in x-direction of the cuboid

dim_y = 0.6 # width of the cuboid

dim_z = 0.7 # height of the cuboid

# initial position of the penetrating sphere

sphere_x = (1/2) * dim_x

sphere_y = 0

sphere_z = dim_z + sphere_r + node_r

engine_gravity = (0,0,0) # gravity acceleration

engine_damping = 0.0 # general damping factor

plot_label = "plotter"

plotData_period = 51 # iterations, uneven for changing nr of interactions

"""" The simulation engine """

O.engines=[

ForceResett

InsertionSo

Interaction

),

NewtonInteg

PyRunner( command = 'addPlotData()', iterPeriod = plotData_period),

]

""" Materials """

O.materials.append(

CohFrictMat( young = node_E, poisson = node_poisson, density = node_rho,

label = node_mat_name

)

)

O.materials.append(

FrictMat(

young = pfacet_E, poisson = pfacet_poisson, density = pfacet_rho,

label = pfacet_mat_name

)

)

O.materials.append(

FrictMat(

young = sphere_E , poisson = sphere_poisson, density = sphere_rho,

label = sphere_mat_name

)

)

""" Bodies """

# adding the 8 corner nodes of the cuboid

nodesIds = [] # holds all gridnode ids

nodesIds.append( O.bodies.append( gridNode( [0,0,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,0,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [0,dim_y,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,dim_y,0], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [0,0,dim_z], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,0,dim_z], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [0,dim_y,dim_z], node_r, wire=False, fixed=True, material=

nodesIds.append( O.bodies.append( gridNode( [dim_x,

# adding the connections and pfacets with the pfacetCreator3 function

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

pfacetCreator3(

# the penetrating sphere

idSphere = O.bodies.append( sphere( [sphere_

""" Movements """

O.bodies[

""" PyRunner Functions """

def addPlotData() :

# if len( O.bodies[

# penetrationDepth = 0

# force_ball = 0

# if O.iter > 200000 :

# O.save(workingDir + '/' + sim_save_folder + filename + '.yade.bz2')

# plot.saveGnuplo

# O.pause()

# else :

# penetrationDepth = O.bodies[

# force_ball = O.bodies[

plot.addData(

time = O.time,

nrInt = O.interactions.

# nrInt = len( O.bodies[

# nrInt = len( O.interactions),

# penetrationDepth = penetrationDepth,

# force_ball = force_ball,

pos_ball = O.bodies[

vel_ball = O.bodies[

forceZ_ball = O.forces.

)

""" Plots """

# plot.plots = {'penetrationDepth' : 'force_ball', 'time' : 'penetrationDepth', 'time ' : 'force_ball', 'time ' : 'nrInt'}

plot.plots = {'time' : 'nrInt', 'time ' : 'pos_ball', 'time ' : 'vel_ball', 'time ' : 'forceZ_ball'}

plot.plot( subPlots = True )

""" Viewer """

qt.Controller()

qtv = qt.View()

qtv.ortho = True

qtv.eyePosition = [0., 1.5, 0]

qtv.viewDir = [0, -1, 0]

qtv.axes = False

qtv.grid = (False, False, False)

qtr = qt.Renderer()

qtr.light2=True

qtr.lightPos=

qtr.bgColor=[1,1,1]

## Can you help with this problem?

Provide an answer of your own, or ask Paul Pircher for more information if necessary.