illogical outputted data for 2D uniaxial compression study

Asked by Kieran

I've adapted an open source Yade code for uniaxial compression of particles in a cube. For reference, the file can be found here: https://gitlab.com/yade-dev/trunk/blob/master/doc/sphinx/tutorial/03-oedometric-test.py

I've altered the code (which is copied below) for the following needs:
* The 3D cube has been reduced to a 2D simulation. This was done by slimming the width of the cube to equal that of the diameter of all of the particles (radius = constant).
* I've replaced the opensource code that tracks the unbalanced force evolution to record the Cartesian coordinates of each particle during compression.

To make an easier trouble shooting process, the original 03-oedometric-test.py code has been commented out and replaced by my choosing.

The code I am running can be seen here:
 # 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
readParamsFromTable(rMean=.03,rRelFuzz=0,maxLoad=1.5e4,minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((0.5,0,0.5),(.5,.03,.5),wallMask=31))
#O.bodies.append(geom.facetBox((0,0,0),(.5,.03,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,0.03,1),rMean=rMean,rRelFuzz=rRelFuzz)
#sp.makeCloud((0,0,0),(-0.5,0.03,-0.25),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation()

O.engines=[
 ForceResetter(),
 # sphere, facet, wall
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
 # the label creates an automatic variable referring to this engine
 # we use it below to change its attributes from the functions called
 PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.5*PWaveTimeStep()

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time

# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
 # 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<5000: return
 # the rest will be run only if unbalanced is < .1 (stabilized packing)
 if unbalancedForce()>.1: return
 # add plate at the position on the top of the packing
 # the maximum finds the z-coordinate of the top of the topmost particle
 O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=0))
 global plate # without this line, the plate variable would only exist inside this function
 plate=O.bodies[-1] # the last particles is the plate
 # Wall objects are "fixed" by default, i.e. not subject to forces
 # prescribing a velocity will therefore make it move at constant velocity (downwards)
 plate.state.vel=(0,0,-.1)
 # start plotting the data now, it was not interesting before
 O.engines=O.engines+[PyRunner(command='myaddData()',iterPeriod=200)]
 # next time, do not call this function anymore, but the next one (unloadPlate) instead
 checker.command='unloadPlate()'

#save parameters of spheres in file

def myaddData():
   if (O.iter%8==0):
      # open file
      fil=open("timeStep_"+str(O.iter)+'.txt',"a+")
      #fil=open("timeStep_"+str(O.iter),"a+")
      #fil.write ('id \t x \t y \t z \t vz \t radius \t force \n')
      for b in O.bodies:
         if (type(b.shape)==Sphere):
            number,rad=b.id,b.shape.radius
            x,y,z=b.state.pos[0],b.state.pos[1],b.state.pos[2]
            veloc,force=b.state.vel[2],O.forces.f(b.id)
            fil.write("{:d}\t{:.3g}\t{:.3g}\t{:.3g}\t{:.3g}\t{:.3g}\t{}\n".format(number,x,y,z,veloc,rad,force)) # ( .format = % )

def unloadPlate():
 # if the force on plate exceeds maximum load, start unloading
 if abs(O.forces.f(plate.id)[2])>maxLoad:
  plate.state.vel*=-1
  # next time, do not call this function anymore, but the next one (stopUnloading) instead
  checker.command='stopUnloading()'

def stopUnloading():
 if abs(O.forces.f(plate.id)[2])<minLoad:
  # O.tags can be used to retrieve unique identifiers of the simulation
  # if running in batch, subsequent simulation would overwrite each other's output files otherwise
  # d (or description) is simulation description (composed of parameter values)
  # while the id is composed of time and process number
  plot.saveDataTxt(O.tags['d.id']+'.txt')
  O.pause()

#def addPlotData():
# if not isinstance(O.bodies[-1].shape,Wall):
# plot.addData(); return
# Fz=O.forces.f(plate.id)[2]
# plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots={'i':('unbalanced',),'w':('Fz',)}
plot.plot()

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()

I have (2) issues with my code:

(1)
The .py file seems to randomly generate the Cartesian coordinates upon simulation with low probability. The norm that I have experienced is the .py file will resort to the original command previously seen in the opensource "03-oedometric-test.py" where it will output empty .txt files with the following information:

# Fz i unbalanced w

However, on random occasions the model will output .txt files that have actually logged the Cartesian coordinates of each particle. This occurs very rarely. I am not sure what causes the program to do one or the other, but it seems I have no real control on how to generate the data.

(2)
I've began data processing on one of the few "successful?" simulations that actually outputted particle Cartesian coordinates. However, the data seems off. To reiterate, the width of the cube in the y-axis is equal to the constant radii of all particles in the simulation. This is how I've simplified a complex 3D study to a 2D simulation for simple data processing. However, some particle y-coordinates are larger than expected with respect to the volume confinement. This makes the data illogical because it is suggesting that some particles are sitting side-by-side in the y-axis, which is limited to the size of one particle.

Any suggestions on how to solve problems (1) and (2)?

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

1)
> PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker')

do not use realPeriod, use iterPeriod or virtPeriod instead.
realPeriod means real Earth time, which is "random" w.r.t. simulation by definition (is different on different machines, on the same machine it is influenced by other processes, RAM usage etc.). This is (probably) the reason why you get "random" output.

cheers
Jan

Can you help with this problem?

Provide an answer of your own, or ask Kieran for more information if necessary.

To post a message you must log in.