Experimental behaviour is different from simulated one

Asked by Luca Strobino on 2018-12-02

Hi everyone,
I've done some tests in laboratory and I'm trying to reproduce the experimental behaviour using YADE but there is something that I don't understand.
The geometry of the problem is the following:
There is an hollow aluminum cylinder fixed to the ground and filled with 20000 aluminum alloy spheres (granular material). An object (named "intruder") , almost in contact with the granular material, is let free to fall inside the cylinder with a given initial velocity ( lets say impact velocity). The objective is to dissipate its kinetic energy by means of the granular material. The problem is that the granular material seems to dissipate the kinetic energy in really more efficient way than what happens in laboratory.
I'm using the following values for the materials:

idIntruder = O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.3,density = 2700,label="intruder"))
idCylinder = O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.3,density = 2700,label="cylinder"))
idSphere = O.materials.append(FrictMat(young=72e9,poisson=.33,frictionAngle=.3,density = 2690,label="sphere"))

and the engine is like this:

O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb(),
                                              Bo1_Box_Aabb(),
                                              Bo1_Polyhedra_Aabb()]),
                           InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),
                                              Ig2_Box_Sphere_ScGeom(),
                                              Ig2_Facet_Polyhedra_PolyhedraGeom(),
                                              Ig2_Wall_Polyhedra_PolyhedraGeom(),
                                              Ig2_Sphere_Polyhedra_ScGeom()],
                                              [Ip2_FrictMat_FrictMat_FrictPhys()],
                                              [Law2_ScGeom_FrictPhys_CundallStrack()]),
                           NewtonIntegrator(gravity=(0,0,-9.81)),
                           VTKRecorder(iterPeriod=37000,recorders=['spheres','boxes','force','stress'],fileName=O.tags['id']+'/'+'vtk'),
                           PyRunner(command="vtkExporter.exportSpheres(ids='all',what=['Sphere_forces','abs(O.forces.f(b.id))'])",iterPeriod=1000),
                           PyRunner(command="vtkExporter.exportPolyhedra(ids='all')",iterPeriod=37000),
                           PyRunner(command='stopper()',iterPeriod=1000),
                           PyRunner(command='addPlotData()',iterPeriod=1000)]

If more information are needed do not esitate to ask please, it's the first time I write here and I'm not sure if I explained the problem in an understandable way. It have been days that I'm trying to understand where the issue could be and maybe someone more expert than me will find it quickly. I have very short time to deliver my work and this forum is my last chance. I hope someone can help me.
Thank you so much in advance,
Luca

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2018-12-03
Last reply:
2018-12-19
Jan Stránský (honzik) said : #1

Hello,

for a meaningful help, we will need a complete example [1]. E.g., you are writing about "alloy spheres", but have polyhedral staff in O.engines..

Concerning dissipation, you can try to set
NewtonIntegrator(...,damping=0.0) # [2]
The value is betweeon 0.0 and 1.0 and controls numerical damping.

[1] https://yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.NewtonIntegrator.damping

Luca Strobino (luca.strobino) said : #2

Hi,
thank you so much for your answer. I'll try to answer to what you said:
- There is polyhedral staff in O.engines because the "intruder" has a particular shape (cylindrical tube with at the end a conical head) that I created by means of polyhedra and it has to interact with the spheres.
- I thought that if I don't insert damping in "NewtonIntegrator" it is set to 0, I'm i right or it's necessary to explicitly set it to 0?
- I'm not sure of what you mean with a complete example so I will copy the code here ( feel free to make all the considerations or to say any suggestions please):

from yade import pack, plot
from numpy import linspace
from yade import ymport, export

#-----------------------------------------------------------------------------#
# MATERIALS DEFINITION #
#-----------------------------------------------------------------------------#

# Materials of the intruder, of the cone and of the spheres

#idAluminium =O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.1,density = 2700,label="aluminium"))
idIntruder = O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.3,density = 2700,label="intruder"))
idCylinder = O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.3,density = 2700,label="cylinder"))
idSphere = O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.3,density = 2700,label="sphere"))
idSatellite = O.materials.append(FrictMat(young=80e9,poisson=.33,frictionAngle=.3,density = 2800,label="satellite"))

#-----------------------------------------------------------------------------#
# TEST MASS DEFINITION #
#-----------------------------------------------------------------------------#

satellite_average_density = 2800.0 #[kg/m^2]
additional_load = 12+(117 - 109.36319) + (117-115.35958797367785) #[kg]
satellite_cube_mass = 94.62250930587025 + additional_load #[kg]
satellite_volume = satellite_cube_mass/satellite_average_density #[m^3]
satellite_radius = ((satellite_cube_mass/satellite_average_density)**(1.0/3.0))/2.0 #[m]

# N.B satellite_radius is half of the latus of the cube that represents the satellite

# Variable to choose the initial heigh (if jump = 0 the initial heigh is the top of the cylinder, with + starts inside with - outside)

granular_level = 0 - 0.0819
jump = granular_level

#-----------------------------------------------------------------------------#
# CYLINDER DEFINITION #
#-----------------------------------------------------------------------------#

cylinder_heigh = 0.5 # [m]
cylinder_radius = 0.0526 # [m]
cylinder_thickness = 0.005 # [m]
cylinder_bottom_thickness = 0.02 # [m]
cylinder_upper_thickness = 0.02 # [m]

granular_heigh = 0.29 # [m]

# Color definition:

intr_color = 0.792
sat_color = 0.6
bottom_color = 0.35
wall_color = 0.88
#-----------------------------------------------------------------------------#
# ELEMENTS CONNECTION #
#-----------------------------------------------------------------------------#

# Here all the elements (satellite + intruder) are clamped together

bodyList = []

#-----------------------------------------------------------------------------#
# INTRUDER CREATION #
#-----------------------------------------------------------------------------#

intruder_radius = 0.028
intruder_heigh = 0.5
intruder_head_heigh = 0.02

number_element_intruder = 30
tethas = linspace(0,2*pi,num = number_element_intruder,endpoint=True)
plus = (2.0*pi/30.0)
rad = intruder_radius

# N.B Reducing or increasing the number "num" in the variable thetas it changes
# also the value of the mass of the intruder. By hand, the value of the mass of
# the intruder is 3.364 [kg]

for th in tethas: bodyList.append(O.bodies.append(utils.polyhedron([(0,0,cylinder_heigh/2 - intruder_head_heigh - jump),
                                                                    (0,0,cylinder_heigh/2 - intruder_head_heigh + intruder_heigh - jump),
                                                                    (rad*cos(th),rad*sin(th),cylinder_heigh/2 -jump),
                                                                    (rad*cos(th + plus),rad*sin(th + plus),cylinder_heigh/2 + intruder_heigh -jump)],
                                                                    color=(intr_color,intr_color,intr_color),
                                                                    material="intruder")))

#-----------------------------------------------------------------------------#
# SATELLITE CREATION #
#-----------------------------------------------------------------------------#

bodyList.append(O.bodies.append(utils.box((0,0,cylinder_thickness + intruder_heigh + intruder_heigh/2 + satellite_radius - jump),
                                      (satellite_radius,satellite_radius,satellite_radius),
                                      color=(sat_color,sat_color,sat_color),
                                      material="satellite")))

#-----------------------------------------------------------------------------#
# UNION OF ALL ELEMENTS #
#-----------------------------------------------------------------------------#

# Piston is formed by test mass and intruder

piston = O.bodies.clump(bodyList)
satellite = O.bodies[piston]

#-----------------------------------------------------------------------------#
# SATELLITE DOF BLOCKING #
#-----------------------------------------------------------------------------#

for i in O.bodies: i.state.blockedDOFs = 'xyXYZ'
for i in O.bodies: i.state.vel=(0,0,-1.14)

#-----------------------------------------------------------------------------#
# CYLINDER CREATION #
#-----------------------------------------------------------------------------#

def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):
   """creates cylinder made of boxes. Axis is parallel with z+"""
   center = Vector3(center)
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles]
   ret = []
   for i,pt in enumerate(pts):
      l = pi*radius/nSegments
      es = (.5*thick,l,.5*height)
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      ret.append(box(pt,es,ori,**kw))
   return ret

number_element_cylinder = 50

surf = cylSurf((0,0,cylinder_bottom_thickness),
               cylinder_radius,
               cylinder_heigh,
               nSegments=number_element_cylinder,
               thick=0.005,
               wire=True,
               material = "cylinder",
               color=(sat_color,sat_color,sat_color))

O.bodies.append(surf)

for b in surf: b.state.blockedDOFs='xyzXYZ'

#-----------------------------------------------------------------------------#
# CYLINDER BOTTOM THICKNESS #
#-----------------------------------------------------------------------------#

cylinder_lower_boundary =O.bodies.append(utils.box((0,0,cylinder_bottom_thickness/2),
                                                   (0.055,0.055,cylinder_bottom_thickness/2),
                                                   material="cylinder",
                                                   color=(bottom_color,bottom_color,bottom_color)))
O.bodies[cylinder_lower_boundary].state.blockedDOFs = 'xyzXYZ'

#-----------------------------------------------------------------------------#
# SPHERE PACK IMPORTATION #
#-----------------------------------------------------------------------------#

packing = ymport.text('gravity_deposition_2_5mm_20000_v3.txt',scale=1.0,material="sphere")
granuli = O.bodies.append(packing)

#-----------------------------------------------------------------------------#
# FLOOR CREATION #
#-----------------------------------------------------------------------------#

ground = utils.wall(0,axis=2,sense=200,color=(0.255,0.255,0))
O.bodies.append(ground)

#-----------------------------------------------------------------------------#
# STROKE DEFINITION #
#-----------------------------------------------------------------------------#

corrective_factor = 0.0027809013496153945

initial_heigh = (cylinder_bottom_thickness + granular_heigh + intruder_head_heigh + intruder_heigh + satellite_radius) + corrective_factor

#-----------------------------------------------------------------------------#
# MASSES DEFINITION #
#-----------------------------------------------------------------------------#

cylinder_mass = b.state.mass * number_element_cylinder
cylinder_lower_boundary_mass = ((0.058*2) * (0.058*2) * cylinder_bottom_thickness) * 2700
granular_material_mass = 0.0001767 * 20000
intruder_mass = satellite.state.mass - satellite_cube_mass

#-----------------------------------------------------------------------------#
# VTK DEFINITION #
#-----------------------------------------------------------------------------#

os.mkdir(O.tags['id'])
vtkExporter = export.VTKExporter(O.tags['id']+'/'+'vtk_files')
vtkExporter.exportPolyhedra(ids = 'all')

#-----------------------------------------------------------------------------#
# ENGINE DEFINITION #
#-----------------------------------------------------------------------------#

# Engine of the simulation

O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb(),
                                              Bo1_Box_Aabb(),
                                              Bo1_Polyhedra_Aabb()]),
                           InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),
                                              Ig2_Box_Sphere_ScGeom(),
                                              Ig2_Facet_Polyhedra_PolyhedraGeom(),
                                              Ig2_Wall_Polyhedra_PolyhedraGeom(),
                                              Ig2_Sphere_Polyhedra_ScGeom()],
                                              [Ip2_FrictMat_FrictMat_FrictPhys()],
                                              [Law2_ScGeom_FrictPhys_CundallStrack()]),
                           NewtonIntegrator(gravity=(0,0,-9.81)),
                           VTKRecorder(iterPeriod=37000,recorders=['spheres','boxes','force','stress'],fileName=O.tags['id']+'/'+'vtk'),
                           #PyRunner(command="vtkExporter.exportSpheres(ids='all',what=['Sphere_forces','abs(O.forces.f(b.id))'])",iterPeriod=1000),
                           #PyRunner(command="vtkExporter.exportFacets(ids='all')",iterPeriod=1000),
                           PyRunner(command="vtkExporter.exportPolyhedra(ids='all')",iterPeriod=37000),
                           #PyRunner(command="vtkExporter.exportBoxes(ids='all')",iterPeriod=1000),
                           PyRunner(command='stopper()',iterPeriod=1000),
                           PyRunner(command='addPlotData()',iterPeriod=1000)]

#-----------------------------------------------------------------------------#
# STOPPER FUNCTION DEFINITION #
#-----------------------------------------------------------------------------#

# With this command I define a stopper function: when the time arrives at n
# seconds the simulation stops and the data are saved in a .txt file

def addPlotData():
 from yade import plot
 plot.addData(time = O.time, satellite_velocity = satellite.state.vel[2], stroke = initial_heigh - satellite.state.pos[2])

def stopper():
 from yade import export
 if O.time > 1.2:
  O.pause()
  plot.saveDataTxt('experiment_result_v3_065_03.txt',vars=('time','satellite_velocity','stroke'))

# This command saves the simulation, so that it can be reloaded later, for experimentation

O.dt = utils.PWaveTimeStep()

O.saveTmp()
O.run()
utils.waitIfBatch()

P.s I hope that everything is clear. Feel free to ask for more information if you need.

Jan Stránský (honzik) said : #3

> I thought that if I don't insert damping in "NewtonIntegrator" it is set to 0, I'm i right or it's necessary to explicitly set it to 0?

no, the default value is 0.2. So please try it with 0.0 and let us know

> I'm not sure of what you mean with a complete example

the code which we can copy-paste and directly test. So apart from the code itself, we would need also "gravity_deposition_2_5mm_20000_v3.txt" file. But before sending it (rather large file I guess), please try the damping.

Jan

Luca Strobino (luca.strobino) said : #4

Thank you so much, I'm running some tests in order to understand if setting damping = 0 there are significant changes. I'll let you know.

Launchpad Janitor (janitor) said : #5

This question was expired because it remained in the 'Open' state without activity for the last 15 days.