Asked by ehsan benabbas on 2020-02-01

Hello everyone,

I want to define a packing of spheres with the following properties:
r_min=0.0001 (m)
r_avr=0.0002 (m)
r_max=0.0003 (m)
number of particles = 20000 (at least, It can be more than this)
cube size = 0.02 * 0.02 * 0.02 (m3)

I am using the makeCloud to do so, the problem is when I print the minimum and maximum value of radii, I realized they have been changed because of "internalCompaction=True".
My goal is to not decrease the sample dimensions, not increase the particles' radii, just increase the number of particles to get the sample and target porosity. Which obviously my code works based on increasing the radius now.

Based on the document of YADE, there are 2 options for packing which are moving the walls or increasing the size of particles that we can implement those by maxWallVelocity and finalMaxMultiplier functions. My desire is to not change these 2 and instead, increase the number of particles. Is that possible to do so?

My code:

######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS #########

print ('************** START **************')
import numpy as np
import time
import datetime, os
start_time=time.time()
from datetime import datetime
import math
from yade import qt, export, utils

######### DEFINING VARIABLES #########

print ('============ DEFINING VARIABLES ============')

num_spheres=20000,
compFricDegree = 29,
key='_triax_',
unknownOk=True
)

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.4
compFricDegree = table.compFricDegree
finalFricDegree = 29
IP=100 # iteration period to record data and stuff
micro_record_iterPeriod=IP
ORN=3000 # O.Run Number of iterations
micro_record_enable_normal_branch=True
micro_record_float_type = np.float32
damp=0.2
thick=0
stabilityThreshold=0.01
PCPC=0.0001 # Precision of Confining Pressure Convergence
r_min=0.1*1e-3 # m
d_min=2*r_min # m
r_max=0.3*1e-3 # m
d_max=2*r_max # m
r_avr=(r_min+r_max)/2 # m
d_avr=2*r_avr # m
r_fuz=(r_max/r_avr)-1 # m
Kn=10e8*(d_avr) ### FIXME
Kt=10e8*(d_avr) ### FIXME
young=Kn/r_avr # 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr
#young=5e6 # contact stiffness
poisson=Kn/Kt # Kt/Kn
#poisson=0.3 # Kt/Kn
Ls=0.02 # m length of specimen ### FIXME
L_REV=7*(d_avr) # m
if Ls < L_REV:
sys.exit("*** ERROR! The specimen's dimension is too samll! ***")
elif Ls==L_REV:
print ("*** This is the minimum specimen's dimension you can take! ***")
else:
print ("*** The specimen's dimension is good enough! ***")
mn,mx=Vector3(0,0,0),Vector3(Ls,Ls,Ls)
Vt=-1*1e-3 # m/s # negative sign describes the compression direction
strainRate=Vt/Ls
#strainRate=-0.01 # %/sec ### FIXME
target_strain=0.2 ### FIXME
print ("The target strain has been set to:", target_strain)
sigmaIso=-5e5 # Pa ### FIXME
particleDensity=2000 #kg/m3
#particleDensity=2600

################# DEFINING MATERIALS #################

print ('============ DEFINING MATERIALS ============')
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))

####################################################
################# DEFINING PACKING #################

print ('============ DEFINING PACKING ============')
walls=aabbWalls([mn,mx],thickness=thick,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
clumps=False
sp.makeCloud(mn,mx,r_avr,r_fuz,num_spheres,False, 0.95,seed=1)
#sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same

os.mkdir('3axresults')
os.chdir('/home/ehsan/Desktop/3axresults')
export.text('InitialPackingData')

################# DEFINING TRIAXIAL TEST #################

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/young,
thickness = thick,
internalCompaction=True,
)

################# DEFINING FUNCTIONS #################

print ('============ DEFINING FUNCTIONS ============')
def history():
e11 = -triax.strain,
e22 = -triax.strain,
e33 = -triax.strain,
ev = -triax.strain-triax.strain-triax.strain,
s11 = -triax.stress(triax.wall_right_id),
s22 = -triax.stress(triax.wall_top_id),
s33 = -triax.stress(triax.wall_front_id),
i = O.iter,
t = O.time, # virtual (yade) time --- time of simulation
fab = utils.fabricTensor())

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

print ('============ DEFINING ENGINES ============')
newton=NewtonIntegrator(damping=damp)

O.engines=[
ForceResetter(),
#GravityEngine(gravity=(0,-9.806,0),warnOnce=False),
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,
PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'),
TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key),
newton
]

Gl1_Sphere.stripes=True

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

print ('============ APPLYING CONFINING PRESSURE ============')

triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = sigmaIso

while 1:
O.run(ORN,True)
unb = unbalancedForce()
meanS=(triax.stress(triax.wall_right_id)+triax.stress(triax.wall_top_id)+triax.stress(triax.wall_front_id))/3
ConfStressRatio=abs(sigmaIso-triax.meanStress)/abs(sigmaIso)
print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
print ('mean stress engine:',triax.meanStress,' mean stress (Calculated):',meanS, ' ConfSratio:',ConfStressRatio,' step:', O.iter/ORN, ' Time:',O.time, ' TimeStep',O.dt)
print ('porosity:',triax.porosity, ' void ratio:',triax.porosity/(1-triax.porosity))
if unb<stabilityThreshold and ConfStressRatio<PCPC:
break

export.text('FinalPhase01PackingData')
e22Check=-triax.strain ###%%%***
print ('Axial Strain',e22Check)
print ('Mean stress engine: ',triax.meanStress)
print ('Mean stress (Calculated): ',meanS)
print ('################## Isotropic phase is finished and saved successfully ##################')

################# REACHING TARGET POROSITY ##################

print ('============ REACHING TARGET POROSITY ============')
import sys
while triax.porosity>targetPorosity:
compFricDegree = 0.95*compFricDegree
print ('\r Friction: ',compFricDegree,' porosity:',triax.porosity,'step= ',O.iter/ORN,' Time:',O.time, ' TimeStep',O.dt)
sys.stdout.flush()
O.run(ORN,True)

print ('################## Target porosity is reached and compacted state saved successfully ##################')

################ PRINT SOME CHECK VARIABLES #################

print ('Number of elements:', len(O.bodies))
print ('Box Volume engine:', triax.boxVolume)
Vt = (mx-mn)*(mx-mn)*(mx-mn) # total volume of the specimen (box)
print ('Box Volume calculated:', Vt)
if Vt == triax.boxVolume:
print ("*** Volume calculation is Correct. ***")
else:
sys.exit("*** ERROR! Volume calculation is WRONG. ***")
Vs=triax.particlesVolume
print('Total volume of particles (Vs):',Vs)
Vv=Vt-Vs
print('Total volume of voids Calculated (Vv):',Vv)
print ('porosity:',triax.porosity)
n=Vv/Vt
print ('porosity Calculated (n):',n)
if n == triax.porosity:
print ("*** Porosity calculation is Correct. ***")
e=n/(1-n)
print ('Void ratio Calculated (e):',e)
else:
sys.exit("*** ERROR! Porosity calculation is WRONG. ***")

Output print:

Friction: 6.552130688784401 porosity: 0.4010231201202055 step= 634.0 Time: 0.13542494035980787 TimeStep 1.1660813485642304e-07
################## Target porosity is reached and compacted state saved successfully ##################
Number of elements: 20006
Box Volume engine: 8.000000000000001e-06
Box Volume calculated: 8.000000000000001e-06
*** Volume calculation is Correct. ***
Total volume of particles (Vs): 4.800954226807098e-06
Total volume of voids Calculated (Vv): 3.1990457731929034e-06
porosity: 0.39988072164911287
porosity Calculated (n): 0.39988072164911287
*** Porosity calculation is Correct. ***
Void ratio Calculated (e): 0.6663354037683562

Thank you,
Ehsan

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
2020-02-02
Last query:
2020-02-02
2020-02-02 Robert Caulk (rcaulk) said on 2020-02-02: #1

>>My goal is to not decrease the sample dimensions, not increase the particles' radii, just increase the number of particles to get the sample and target porosity.

If cube size cannot change, and the distribution of r cannot change, then the only thing you can do is add more spheres during the make cloud step . I hope it is obvious that you cannot apply stresses without changing particle size or box dimensions.

 ehsan benabbas (ehsanben) said on 2020-02-02: #2

Thanks Robert Caulk, that solved my question.