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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp