yade-users team mailing list archive

[Paul Pircher <question696480@xxxxxxxxxxxxxxxxxxxxx>]

New question #696480 on Yade:
https://answers.launchpad.net/yade/+question/696480

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=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Node_Aabb(), Bo1_GridConnection_Aabb(), Bo1_PFacet_Aabb(), ]),
    InteractionLoop(
        [Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
        Ig2_Sphere_PFacet_ScGridCoGeom(), Ig2_PFacet_PFacet_ScGeom(),
        Ig2_Sphere_GridConnection_ScGridCoGeom(), Ig2_Sphere_Sphere_ScGeom(), ],           
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
        Ip2_FrictMat_FrictMat_FrictPhys() ],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), Law2_ScGeom_FrictPhys_CundallStrack(), 
        Law2_ScGridCoGeom_FrictPhys_CundallStrack(), Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), ]
    ),
    NewtonIntegrator(  gravity=engine_gravity,  damping=engine_damping),
    PyRunner( command = 'addPlotData()', iterPeriod = plotData_period),
]

""" Materials """
O.materials.append( 
    CohFrictMat( young = node_E, poisson = node_poisson, density = node_rho,
        frictionAngle = radians(node_phi), normalCohesion = node_normal_coh,
        shearCohesion = node_shear_coh, momentRotationLaw = True,
        label = node_mat_name
    ) 
)
O.materials.append( 
    FrictMat( 
        young = pfacet_E, poisson = pfacet_poisson, density = pfacet_rho,
        frictionAngle = radians(pfacet_phi),
        label = pfacet_mat_name 
    )
) 
O.materials.append( 
    FrictMat( 
        young = sphere_E , poisson = sphere_poisson, density = sphere_rho,
        frictionAngle = radians(sphere_phi),
        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='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [dim_x,0,0], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [0,dim_y,0], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [dim_x,dim_y,0], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [0,0,dim_z], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [dim_x,0,dim_z], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [0,dim_y,dim_z], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) )
nodesIds.append( O.bodies.append( gridNode( [dim_x,dim_y,dim_z], node_r, wire=False, fixed=True, material='gridNodeMat' ) ) ) 

# adding the connections and pfacets with the pfacetCreator3 function
pfacetCreator3(0,1,3, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(0,2,3, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(4,5,6, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(5,6,7, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(0,1,4, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(1,4,5, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(2,3,7, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(2,6,7, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(1,3,7, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(1,5,7, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(0,2,4, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)
pfacetCreator3(2,4,6, cylIds=[], pfIds=[], wire=False, material=pfacet_mat_name, color=color)

# the penetrating sphere 
idSphere = O.bodies.append( sphere( [sphere_x,sphere_y,sphere_z] ,sphere_r,wire=False,fixed=False,material = sphere_mat_name) )

""" Movements """
O.bodies[idSphere].state.vel = (0,0, sphere_v) # initial velocity of the sphere

""" PyRunner Functions """
def addPlotData() :
    # if len( O.bodies[idSphere].intrs() ) == 0 :
    #     penetrationDepth = 0
    #     force_ball = 0
        # if O.iter > 200000 :
        #     O.save(workingDir + '/' + sim_save_folder + filename + '.yade.bz2')
        #     plot.saveGnuplot(data_save_folder + filename) 
        #     O.pause()
    # else :
    #     penetrationDepth = O.bodies[idSphere].intrs()[0].geom.penetrationDepth * 1e3
    #     force_ball = O.bodies[idSphere].intrs()[0].phys.normalForce.norm()
    plot.addData(
        time = O.time,
        nrInt = len( O.bodies[idSphere].intrs() ),
        # penetrationDepth = penetrationDepth,
        # force_ball = force_ball,
        pos_ball = O.bodies[idSphere].state.pos[2],
        vel_ball = O.bodies[idSphere].state.vel[2],
        forceZ_ball = O.forces.f(idSphere)[2],
    )

""" 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=Vector3(1200,1500,500)
qtr.bgColor=[1,1,1]

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.