Target porosity in spherepack.makecloud

Asked by Shyam Kasi

Hi, I was trying to do a triaxial test on the sand of porosity 20% and I wasn't able to reach the input porosity. I used the experimental data to determine the particle size distribution and here below I am presenting the data.
Sieve Size(mm) Weight retained on the sieve(grams)
0 0
0.075 1.94
0.15 2.26
0.212 155.33
0.425 303.57
0.6 35.3
The total weight of the sand is 500 grams.
I tried changing the distribution, other solutions like reducing friction angle method none of those achieved the required porosity. I am always getting a porosity of more than 40% and reduced to 37% by using the friction method.

Can anyone tell me if I entered the psdSize and psdComm array correctly in the below line of code?

sp.makeCloud(mn,mx,psdSizes=[0,0.000075,0.00015,0.000212,0.000425,0.0006],psdCumm=[0,3.88e-3,8.4e-3,0.3190,0.9262,1],porosity=0.2)

Question information

Language:
English Edit question
Status:
Needs information
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#1

Hi Kasi,

next time please provide MWE (minimum working example) [1]. You provided just one line of code, so no one can reproduce your problem without additional assumptions.

However, I can guess that the problem is that you do not use the option 'distributeMass' [2]. Because of that, you do not distribute the mass of the particles but the number. Thus your simulation is dominated by the bigger particles, and it is likely the reason for the problem with achieving the target porosity.

Best wishes,
Karol

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.pack.html?highlight=makecloud#yade._packSpheres.SpherePack.makeCloud

Revision history for this message
Jérôme Duriez (jduriez) said :
#2

Note also that your numerical sample is made of spheres, unlike (I'm assuming) your real reference sand => it's not sure in my opinion you will be able to form packings with the same initial porosity that the experimental one, even though you have the same distribution of particles "size"

Revision history for this message
Robert Caulk (rcaulk) said :
#3

Yes, to add to Jérôme: if your sand is even slightly silty, you will need to account for that first part of the PSD curve more precisely in order to match the real particle size distribution. Problem is that the smallest particles control the stable timestep in your simulation; the small particles prevent any useful simulation in that case.

Revision history for this message
Shyam Kasi (skc18) said :
#4

Thanks for your suggestions.
When I used distributeMass I was getting an error message saying the script is not responding.
I tried changing the particle size distribution but none of those works.
Here below I am presenting my code.

from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=5e6,# number of spheres
 compFricDegree =43.22, # contact friction during the confining phase
 key='_triax_base_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity =0.199 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree =43.22 # contact friction during the deviatoric loading
rate=-0.01 # loading rate (strain rate)
damp=0.4 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=58e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(0.005,0.005,0.005) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2655.6,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

clumps=False #turn this true for the same example with clumps
if clumps:
 ## approximate mean rad of the futur dense packing for latter use
 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 ## define a unique clump type (we could have many, see clumpCloud documentation)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 ## generate positions and input them in the simulation
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()#get more accurate clump masses/volumes/inertia
else:
 sp.makeCloud(mn,mx,psdSizes=[0,0.000075,0.00015,0.000212,0.000425,0.0006],psdCumm=[0,3.88e-3,403.88e-3,0.7184,0.95554,1],porosity=0.2,distributeMass=True) #"seed" make the "random" generation always #the same
 #sp.makeCloud(mn,mx,psdSizes=[0,0.000075,0.00015,0.000212,0.000425,0.0006],psdCumm=[0,0.1,0.3,0.62,0.92554,1],porosity=0.2)
 #sp.makeCloud(mn,mx,rMean=0.00333/2,rRelFuzz=0,num=2000,porosity=0.2)
 #O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
 #or alternatively (higher level function doing exactly the same):
 sp.toSimulation(material='spheres')

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 #maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 #finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=False,
 max_vel=0.001 # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(gravity=(0,-9.81,0),damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]

#Display spheres with 2 colors for seeing rotations better
#Gl1_Sphere.stripes=0
#if nRead==0: yade.qt.Controller(), yade.qt.View()

## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE
## DEPENDING ON YOUR EDITOR, IT COULD BE DONE
## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES
## AND PRESSING CTRL+SHIFT+D

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-50000

while 1:
 O.run(1000, True)
   ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
   unb=unbalancedForce()
   print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
   if unb<stabilityThreshold and abs(-50000-triax.meanStress)/50000<0.01:
      break

O.save('confinedState'+key+'.yade.gz')
print "### Isotropic state saved ###"
O.pause()
###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

### We will reach a prescribed value of porosity with the REFD algorithm
### (see http://dx.doi.org/10.2516/ogst/2012032 and
### http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
 ## we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 ## while we run steps, triax will tend to grow particles as the packing
 ## keeps shrinking as a consequence of decreasing friction. Consequently
 ## porosity will decrease
 O.run(500,1)

O.save('compactedState'+key+'.yade.gz')
print "### Compacted state saved ###"
O.pause()
##############################
### DEVIATORIC LOADING ###
##############################

##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-50000
triax.goal3=-50000

##we can change damping here. What is the effect in your opinion?
newton.damping=0.4

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
#O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

from yade import plot

### a function saving variables
def history():
 plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
   ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
   i=O.iter)
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=2000,command='history()',label='recorder')]+O.engines[5:7]
#if 1:
     ## include a periodic engine calling that function in the simulation loop
 #O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

 #O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
#else:
   ## With the line above, we are recording some variables twice. We could in fact replace the previous
   ## TriaxialRecorder
   ## by our periodic engine. Uncomment the following line:
 #O.engines[4]=PyRunner(iterPeriod=2000,command='history()',label='recorder')

#O.run(100,True)

### declare what is to plot. "None" is for separating y and y2 axis
plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
### the traditional triaxial curves would be more like this:
plot.plots={'e22':'s22'}

## display on the screen (doesn't work on VMware image it seems)
plot.plot()

##### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) #####

## In that case we can still save the data to a text file at the the end of the simulation, with:
#plot.saveDataTxt('results'+key)
##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot:
#plot.saveGnuplot('plotScript'+key)

Revision history for this message
Jan Stránský (honzik) said :
#5

> I was getting an error

please always provide the complete error message [3]

> error message saying the script is not responding

O.run(N,True) ... True here means wait=True [4], i.e. the script waits until it ends and is not responding because of this waiting.
Most likely it is not an error, but actually desired behavior.

I suppose this "error" is related to GUI window. Is that right?
If so, do not take care about GUI windows.

Cheers
Jan

[3] https://www.yade-dem.org/wiki/Howtoask
[4] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Omega.run

Revision history for this message
Shyam Kasi (skc18) said :
#6

Yeah, that is a pop-up error given by ubuntu. Even after clicking wait, it is taking a very long time to start the confining phase.

Revision history for this message
Robert Caulk (rcaulk) said :
#7

>Even after clicking wait, it is taking a very long time to start the confining phase.

How many particles do you have? What is the maximum and minimum size of the particles?

Revision history for this message
Shyam Kasi (skc18) said :
#8

Around 3000 particles, Max particle diameter is 0.6mm and minimum particle diameter is 0.075mm.

Revision history for this message
Shyam Kasi (skc18) said :
#9

Still facing the same issue even after using "distributeMass=True"

Revision history for this message
Jan Stránský (honzik) said :
#10

"the same issue" means "it is taking a very long time" or something different?

Revision history for this message
Shyam Kasi (skc18) said :
#11

No, now the distributeMass=True is working I mean it is not taking time but the porosity was not achieved.

Can you help with this problem?

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

To post a message you must log in.