# 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:
- Expired

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Last query:

- Last reply: