New question #695154 on Yade: https://answers.launchpad.net/yade/+question/695154
Dear Yade Users, During the simulation of polyhedra deposition in the box, I encountered instability problems (even for small increments, some grains behave like popcorn). I have made some animations to observe this and came up with the conclusion that it may be related to switching the contact conditions. I have read the paper of Eliáš (which is full of clever ideas by the way), and I think that the issue may be the calculation of the direction of the normal force. I conducted a simulation (see the script at the bottom of the post) that seems to confirm that hypothesis. I have two polyhedra, a cube and a prism (half of the cube) that are shifted. I slide one polyhedron to make a contact starting from one side of the cube. As I move the polyhedron, the intersection becomes symmetrical, and later it is "closer" to the other (perpendicular) wall. I would expect the forces to change gradually, but there is a sudden switch between vertical and horizontal forces. Everything below is speculation since I am not a C++ developer. I have found the part of the code responsible for normal force calculation (function FindNormal in [2]). It utilizes a function, that has a known flaw [3], so it is a possible suspect for me. I wish, I could be able to verify this by myself, but it is not my league (yet ;) ). I am a more Python person. [1] Eliáš, J. (2014). Simulation of railway ballast using crushable polyhedral particles. Powder Technology, 264, 458-465. [2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/Polyhedra_support.cpp [3] https://github.com/CGAL/cgal/issues/2676 --------------------------------------------- from yade import polyhedra_utils from yade import plot from yade import qt ######################################## CONSTANTS d = 0.1 ######################################## FUNCTIONS def addPlotData(): f0 = O.forces.f(0)[0] f1 = O.forces.f(0)[1] f2 = O.forces.f(0)[2] t0 = O.forces.t(0)[0] t1 = O.forces.t(0)[1] t2 = O.forces.t(0)[2] plot.addData(t=O.time,f0=f0,f1=f1,f2=f2,t0=t0,t1=t1,t2=t2) ########################## MATERIALS polyMat = PolyhedraMat() ########################## BODIES p1 = polyhedra_utils.polyhedra(polyMat, v= [(0,0,0),(0,d,0),(d,d,0),(d,0,0),(0,0,d),(0,d,d),(d,d,d),(d,0,d)], fixed = True) p1.state.pos = (0,0,0) p2 = polyhedra_utils.polyhedra(polyMat, v= [(0,0,0),(0,d,0),(d/2,d,0),(d/2,0,0),(0,0,d),(0,d,d),(d/2,d,d),(d/2,0,d)], fixed = True) p2.state.pos = (0,0.9*d,d) O.bodies.append((p1,p2)) ########################## ENGINES O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Polyhedra_Aabb()]), InteractionLoop( [Ig2_Polyhedra_Polyhedra_PolyhedraGeom()], [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] ), NewtonIntegrator(), PyRunner(command="addPlotData()",iterPeriod=100), ] O.dt= 0.001*polyhedra_utils.PWaveTimeStep() plot.plots={'t':(('f0','r'),('f1','g'),('f2','b'),('t0','r--'),('t1','g--'),('t2','b--')),} plot.plot(subPlots =False) try: from yade import qt qt.Controller() v = qt.View() v.ortho = True v.viewDir = (-1,0,0) v.showEntireScene() except: pass p2.state.vel = (0,0,-1) O.run(20000) -- 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
