# I can not get the desired porosity with my PSD

Hello All,
I am trying to generate sphere packing according to a given PSD.

I want to get a porosity=0.375 or e=(triax.porosity)/(1-triax.porosity)=0.6

I use num_spheres=1000 and my answer for porosity is 0.4969
I use num_spheres=5000 and my answer for porosity is 0.435
I use num_spheres=15000 and my answer for porosity is 0.411
I use num_spheres=20000 and my answer for porosity is 0.4139
I use num_spheres=30000 and my answer for porosity is 0.4059
I use num_spheres=40000 and my answer for porosity is 0.39
I use num_spheres=50000 and my answer for porosity is 0.39
I use num_spheres=60000 and my answer for porosity is 0.39

How to use this psd to get the desired porosity=0.375 or e=(triax.porosity)/(1-triax.porosity)=0.6 ?

=============================================================
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab
from numpy import *
import matplotlib.pyplot as plt
import math

seed=table.seed
num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-2820

## creat a packing with a specific particle side distribution (PSD)
psdSizes,psdCumm=[0.0011,0.0012,0.0030,0.0059,0.0081,0.0111,0.0133,0.0222,0.0285,0.0371,0.0746,0.1491,0.2982,0.4224],[0.,0.22,0.29,0.36,.45,0.55,0.59,0.69,0.76,0.87,0.94,0.95,0.97,1.0]
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(2.5,2.5,2.5)
sp.makeCloud(minCorner=mn,maxCorner=mx,num=num_spheres,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=False,seed=seed)
sp.psd(bins=5000,mass=False)
sp.toSimulation

## create material #0, which will be used as default
O.materials.append(FrictMat(young=15e6,poisson=.4,frictionAngle=0,density=0,label='frictionless'))

## create walls around the packing
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

triax=TriaxialStressController(
internalCompaction=False,
goal1=confiningS,
goal2=confiningS,
goal3=confiningS,
max_vel=10,
label="triax"
)

newton=NewtonIntegrator(damping=0.4)

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()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton
]

O.dt=.1*PWaveTimeStep()
O.dynDt=False

while 1:
O.run(1000, True)
unb=unbalancedForce()
print 'unb: ',unb,' t.p: ',triax.porosity,' n: ',porosity(),' e: ',(triax.porosity)/(1-triax.porosity),' Sigma: ',-triax.stress(3)[1]
if unb<0.001 and abs(confiningS-triax.meanStress)/(-confiningS)<0.001:
break

print 'porosity: ',porosity()
print 'triax.porosity',triax.porosity
print 'e: ',(triax.porosity)/(1-triax.porosity)
==========================================================

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2019-06-06
2019-06-06
 Robert Caulk (rcaulk) said on 2019-06-03: #1

Have you tried changing the PSD to increase the fraction of smaller particles? Have you tried not passing numSpheres to makeCloud? [1]

 Saeed (saeed0668191) said on 2019-06-06: #2

hi Robert
I increase the number of compFricDegree, but porosity (n) and Void Ratio: e=n/(1-n) increased and did not help.
Without using the num_spheres in sp.makeCloud , it simulates the small interval of the psd, and the entire interval of the seeds does not coat, and a very large amount of grain is modeled, which is very long runtime.
And the psd does not fully model the numbers given in psdSizes and psdCumm.
thanks

 Robert Caulk (rcaulk) said on 2019-06-06: #3

Have you tried changing the PSD to increase the fraction of smaller particles?

Have you tried following the triax-example porosity method [1]?

 Bruno Chareyre (bruno-chareyre) said on 2019-06-06: #4

distributeMass=True would probably help