yade-users team mailing list archive

[Luc OGER <question668315@xxxxxxxxxxxxxxxxxxxxx>]

Question #668315 on Yade changed:
https://answers.launchpad.net/yade/+question/668315

    Status: Needs information => Open

Luc OGER gave more information on the question:
Hello Robert
thanks for trying to help me

in more detail, I am filling abox by classical gravity depositon than
slowly titlt the box and look at the sphere displacements.

as the first 'step' made a correct converging force balance, the change of the gravity direction is not playing immediately enough change in the balancedforcecalculation to reach the 'new equilibrium', so I have used the same idea as in your oedometrixc sample case
i hope that I am enough clear?

 So here one lineof the table for my batch run then the full code py. 
table:
pas_box theta_max converg_min nb_layer init_seed
.050 35. .01 20 10

py code:
# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch

# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
# gravity deposition in box, showing how to plot and save history of data,
# and how to control the simulation while it is running by calling
# python functions from within the simulation loop
#pas_box theta_max converg_min
#.05 35. .01
readParamsFromTable(pas_box=0.05, theta_max=35.0, converg_min=0.01,nb_layer=20.0,init_seed=10)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *


# import yade modules that we will use below
from yade import pack, plot, export
  
  
box_size = 5.0
# create rectangular box from facets
O.bodies.append(geom.facetBox((box_size/2.0,box_size/2.0,box_size*0.62),(box_size/2.0,box_size/2.0,box_size*0.62),wallMask=31))

# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
packing_fraction = 60./100.
global  i_pas, pas_box,step0,gravity_y,gravity_z,theta_max,converg_min,nb_layer,init_seed
i_pas = 0 
pas_box=.0250
step0 = 0
nb_layer=20.0
gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)
                    
rayon = 0.025
height = 5. * rayon*nb_layer
nombre= int(packing_fraction*(box_size*box_size*nb_layer)/(rayon*rayon*3.14159) )
print(nombre)
print(i_pas*pas_box)
sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0,0,0),(box_size,box_size,height),rMean=rayon,rRelFuzz=.005,num= nombre,seed=init_seed)
# add the sphere pack to the simulation
sp.toSimulation()


# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 global i_pas, pas_box, step0,gravity_y,gravity_z,theta_max,converg_min
   # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 if O.iter<(step0+5000): return 
 if utils.unbalancedForce()<converg_min:
        i_pas = i_pas+ 1
        gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
        gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)
        Newton_integrator.gravity=Vector3(0,-gravity_y,-gravity_z)
        step0 = O.iter
#        print('%8.3f %8.3f %8.3f %8.3f' step0,pas_box*i_pas,sin( pas_box*i_pas*3.14/180.0),cos( pas_box*i_pas*3.14/180.0))
        print (step0, (pas_box*i_pas), gravity_y, gravity_z)
        if (pas_box*i_pas)>theta_max :
            O.pause()
#        plot.saveDataTxt('bbb.txt.bz2')
        nom_file=str(i_pas)+"_"+str(pas_box)+"_"+str(nb_layer)+"_"+str(converg_min)
        filename1='./output_'+nom_file+'.txt'
        # sauvegarde fichier txt
        export.textExt(filename1,format='x_y_z_r',comment='x y z r')
        filename2='./output_'+nom_file+'.vtk'
        export.text2vtk(filename1,filename2)
        
        
# simulation loop 
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	InteractionLoop(
		# handle sphere+sphere and facet+sphere collisions
	    [Ig2_Sphere_Sphere_ScGeom(),
	     Ig2_Facet_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=radians(15.0))],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),

	NewtonIntegrator(damping=0.4,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
	# call the checkUnbalanced function (defined below) every 2 seconds
	PyRunner(command='checkUnbalanced()',realPeriod=5),
	# call the addPlotData function every 200 steps
	PyRunner(command='addPlotData()',iterPeriod=200)
]
    
#without friction 
O.materials[0].frictionAngle=radians(15.0)
# timesteps
O.dt=.5*PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy=True

    
# collect history of data which will be plotted
def addPlotData():
 global  i_pas, pas_box
# print( pas_box*i_pas,sin( pas_box*i_pas*3.14/180.0),cos( pas_box*i_pas*3.14/180.0)),
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter,unbalanced=utils.unbalancedForce(),**O.energy)

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis
plot.plots={'i':('unbalanced',None,O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()

# save tmp pour rerun convergence
O.saveTmp()

O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.