Experimental behaviour is different from simulated one

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

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.

> 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

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.

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