yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #18491
Re: [Question #676528]: Experimental behaviour is different from simulated one
Question #676528 on Yade changed:
https://answers.launchpad.net/yade/+question/676528
Status: Answered => Open
Luca Strobino is still having a problem:
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.
--
You received this question notification because your team yade-users is
an answer contact for Yade.
