# yade-users team mailing list archive

## [Question #695154]: Polyhedra contact instability (potential Bug?)

 Thread Previous • Date Previous • Date Next • Thread Next

```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.

```

 Thread Previous • Date Previous • Date Next • Thread Next