yade-users team mailing list archive

[Karol Brzezinski <question695154@xxxxxxxxxxxxxxxxxxxxx>]

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.