Unstable behaviour for different Young's modulus'

Hi there,

I've been trying to simulate the behaviour of toppling dominoes, but I've been encountering many problems such as unstable behaviour under certain circumstances. I am a beginner in using YADE and I hope someone can help me.

- So when I want to simulate the 'clumps' with a E-modulus of 8000MPa (so 8e9 Pa in SI), the simulation takes forever, it basically does not topple.
- In order to solve this, I had to use a REALLY low E-modulus (like 800) in order to let it run smoothly and get 'reasonable' results.
                 - However for larger ratio's (so distances), the system randomly collapses due to this very low E modulus.

So my question is, is there a way such that I can use a realistic E-modulus, while having a 'smooth' simulation, which does not take too much time to simulate?

This is my script:

#!/usr/bin/python
# -*- coding: utf-8 -*-

'''Example of modeling dominos in YADE'''
from __future__ import print_function

from builtins import range
from yade import pack,export,qt,plot
from pprint import pprint
import numpy

#define material for all bodies(SI units):
id_Mat=O.materials.append(FrictMat(young=8e2,poisson=0.5,density=532,frictionAngle=0.4))
Mat=O.materials[id_Mat] #Soft wood

id_boxMat=O.materials.append(FrictMat(young=8e9,poisson=0.5,frictionAngle=0.4)) #Some material with relatively high friction to prevent slip

#-----------------------------

# Parameters of a domino

x_step = 0.0035 #Particle size in mm 7*20*44
y_step = 0.005
z_step = 0.0044

N_x = 2 #Particle Distribution
N_y = 4
N_z = 10

dt = x_step * N_x # domino thickness
dw = y_step * N_y # domino width
dh = z_step * N_z # domino height

step = [x_step, y_step, z_step]
N = [N_x, N_y, N_z]

rad = 0.0025 #Radius of a spherical element - should be larger than 1/2 * max[x_step, y_step, z_step] in SI

#create a scene:

num_dominos = 30 #Number of domino bricks
ratio = 0.7 #s/h ratio
spacing = ratio*dh + dt # Spacing between them + thickness of a domino: left point-left point

# Box sizes , The blue surface
ll = 1.5*(num_dominos * spacing) # Lenght
ww = dw + 1*dw #Some extra, arbitrary width given
hh = 0.01

id_box = O.bodies.append(box((ll/3,ww/2.,-hh/2.),(ll/2.,ww/2,hh/2.), fixed=True, color = (0.1,0.4,0.9), material=id_boxMat))

for m in range(num_dominos):
    x_disp = (m+0.5) * spacing
    y_disp = ww/2-dw/2+rad/2
    z_disp = rad

    clump1=O.bodies.appendClumped([sphere([x_disp,0+y_disp,0+z_disp], material=Mat, color = (0.9,0.6,0.9), radius=rad)])
    id_clump1 = clump1[0] #Initial Clump is made here

    l=False
    for i in range(N[0]):
        for j in range(N[1]):
            for k in range(N[2]): #The range of clumps is defined here
                if (l):
                    id_new=O.bodies.append(sphere([i*step[0]+x_disp,j*step[1]+y_disp,k*step[2]+z_disp], material=Mat, color = (0.7,0.9,0.9), radius=rad))
                    O.bodies.addToClump([id_new],id_clump1)
                l = True
#The other particles are defined in order to get the clump

# "Hummer" pushing the 1st domino
frad = 0.008 #Radius of arbitrary hummer in
id_new=O.bodies.append(sphere([spacing/2.-frad-rad,ww/2,dh], material=Mat, fixed = False, color = (0.9,0.1,0.1), radius=frad))

O.dt=1*PWaveTimeStep() # Time integration timestep

O.bodies[-1].state.vel = (0.3,0,0) # Start moving the hummer

# Visualization settings
rr = qt.Renderer()
rr.bgColor = Vector3(1., 1., 1.)
rr.light1 = True
rr.light2 = True
rr.ghosts = False
rr.light2Color = Vector3(0.2,0., 0.)

#define engines:
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], # box-box interactions do not exist in yade, so work with clumps
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(damping=0.22,gravity=[0,0,-9.81]) # Yade damping is rather artificial thing, and everything depends on it!

]
O.engines+=[PyRunner(command='addPlotData()',iterPeriod=10)] #To track data over period
plot.plots={'t ':('x_vel')}
plot.title={'Ratio of %ratio , $\mu$ = %frictionAngle'}
plot.labels={'t':'Virtual time in $s$' , 'x_vel':'Velocity x-direction in m/s'}
def addPlotData():
 sph= id_new
 x_vel=numpy.average([b.state.vel[0] for b in O.bodies])
 plot.addData(t=O.time,i=O.iter,x_vel=x_vel)

plot.plot(subPlots=False)
#plot.liveInterval=.2
#max_vel= max(plot.data['x_vel']) call for max approached velocity

qt.View()
O.run(-1)

Thank you very much,

Daniel