Unstable behaviour for different Young's modulus'

Asked by daniel ding on 2021-04-29

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Karol Brzezinski
Solved:
Last query:
Last reply:
Best Karol Brzezinski (kbrzezinski) said : #1

Hi Daniel,

Probably it is not an answer that you hoped for, but my best answer is to try stiffness somewhere between 8e2 and 8e9. I turned off the movement of your hammer and domino collapsed after 2-2.5 s (with young = 800). However, when I increased stiffness to 8000, the simulation remained stable over 15 virtual seconds (and then I just turned it off). It is more than enough for your simulation (which takes 3 virtual seconds).

Best wishes,
Karol

daniel ding (danielding21) said : #2

Thanks Karol Brzezinski, that solved my question.