Make Radius Ratio for the sample and verity of particles radii

Hello,

I made some changes to Bruno's Triaxial code as follows. The problem is, I want to have the ratio of radii as R_max/R_min=3, however, when I get the output of radii at the end of my code, it seems all particles have the same size as 0.6 mm (as the R_avr is even 0.6 mm as well). Is there anything wrong with my code regarding radius defining or I must use another way to make the sample?

Thank you for your help.

######################################################################################################
######### 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
import shutil
start_time=time.time()
from datetime import datetime
import math
from yade import qt, export, utils
from yade import pack
import yade.timing;
from pprint import pprint
O.timingEnabled=True

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

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

nRead=readParamsFromTable(
 num_spheres=5000,
 compFricDegree = 10,
 key='_triax_',
 unknownOk=True
)

from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.404
compFricDegree = table.compFricDegree
finalFricDegree = 30
IP=2000 # iteration period to record data and stuff
IPmacro=100
micro_record_iterPeriod=IP
ORN=1000 # 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
Ls=0.02 # m
Lh=1*Ls
r_min=0.15e-3 # m
d_min=2*r_min # m
r_max=0.45e-3 # m
d_max=2*r_max # m
r_avr=(r_min+r_max)/2 # m
d_avr=2*r_avr # m
mn,mx=Vector3(0,0,0),Vector3(Ls,Lh,Ls)
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
r_fuz=0.2*r_avr # m
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! ***")
young=2*(1e8) # Kn = 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr
print ('The Young modulus is:', young)
Kn=young*(r_avr) ### FIXME
print ('The Stiffness is:', Kn)
Kt=Kn ### FIXME
poisson=Kt/Kn # Kt/Kn
cvs=1-targetPorosity # coefficient of solid volume ===>>> 1 - n = 1 - (Vv / Vt) = Vs / Vt , cvs*volume=Vs
mean_rad = pow(0.24*cvs*volume/num_spheres,0.3333)
print('mean_rad',mean_rad)
strainRate=-0.5 # %/sec ### FIXME
target_strain=0.4 ### FIXME %
sigmaIso=-5e5 # Pa=N/m2 ### FIXME
particleDensity=2600

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

print ('============ DEFINING FUNCTIONS ============')
from yade import plot
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,
        t = O.time, # virtual (yade) time --- time of simulation
        fab = utils.fabricTensor()[0])

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

print ('============ DEFINING MATERIALS ============')

O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres'))
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')
for w in walls:w.shape.radius=0
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)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

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

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

####################################################
################# 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()]),
    #InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06),
 InteractionLoop(
        #ompThreads=8,
  [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=IPmacro,command='history()',label='macro_recorder'),
    PyRunner(iterPeriod=micro_record_iterPeriod,command='record_micro_data()',label='micro_recorder',dead=True),
 TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key),
 newton
]

###############################################################
################# 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)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
    ConfStressRatio=abs(sigmaIso-triax.meanStress)/abs(sigmaIso)
    print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
    print ('ConfSratio:',ConfStressRatio,' step:', O.iter/ORN, ' Time:',O.time, ' iter/sec',O.speed)
    print ('mean stress engine:',triax.meanStress,' mean stress (Calculated):',meanS)
    if unb<stabilityThreshold and ConfStressRatio<PCPC:
        break

export.text('FinalPhase01PackingData')
e22Check=-triax.strain[1]
print ('Axial Strain',e22Check)
print ('porosity:',triax.porosity)
print ('void ratio:',triax.porosity/(1-triax.porosity))
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:
    O.run(ORN,True)
    compFricDegree = 0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print ('\r Friction: ',compFricDegree,' porosity:',triax.porosity,'step= ',O.iter/ORN,' Time:',O.time, ' iter/sec',O.speed)
    sys.stdout.flush()

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

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

RRmax=max([b.shape.radius for b in O.bodies])
RRmin=min([b.shape.radius for b in O.bodies])
print('R_max:',RRmax)
print('R_min:',RRmin)
print ('Number of elements:', len(O.bodies))
Vt = triax.boxVolume
print ('Vt:', Vt)
Vs = triax.particlesVolume
print('Vs:',Vs)
Vv=Vt-Vs
print('Vv:',Vv)
n = triax.porosity
print ('n:',triax.porosity)
e=n/(1-n)
print ('e:',e)
print ('Step that starts the deviatoric loading ', O.iter/ORN)

######################################################
################# DEVIATORIC LOADING #################

print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))

# Allow walls to move. >>> False means fixed , True means free to move
triax.wall_top_activated=True
triax.wall_bottom_activated=False
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True

# Make Stress Control on X,Y,Z directions by setting True, False means Strain Control
triax.stressControl_1=True
triax.stressControl_2=False
triax.stressControl_3=True

# set the raite of strain to the desirable axis
#triax.strainRate1=100*rate
triax.strainRate2=strainRate
#triax.strainRate3=100*rate

# set the deviatoric conditions
triax.stressMask = 5
triax.goal1 = sigmaIso
triax.goal2 = strainRate
triax.goal3 = sigmaIso

newton.damping=0.1

while 1:
    O.run(ORN,True)
    unb=unbalancedForce()
    axialS=triax.stress(triax.wall_top_id)[1]
    eps2=triax.strain[1]
    eps2cal=(triax.height-Ls)/Ls
    print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
    print ('sigma2: ',axialS, ' q: ', axialS-sigmaIso,' step= ', O.iter/ORN,' Time:',O.time, ' iter/sec',O.speed)
    print ('Strain2 Calculated: ',eps2cal, ' Strain2 engine: ',eps2,'axial deformation (%)', (eps2-e22Check)*100)
    if abs(eps2-e22Check)>=target_strain: #triax.sigma2
        break

print ('################## Deviatoric phase is finished and saved successfully ##################')

#############################################################
#################### SOME DATA TO SEE ######################

RRmax=max([b.shape.radius for b in O.bodies])
RRmin=min([b.shape.radius for b in O.bodies])
RRavr=sum([b.shape.radius for b in O.bodies])/len(O.bodies)
print([b.shape.radius for b in O.bodies])
print('R_max:',RRmax)
print('R_min:',RRmin)
print('R_avr:',RRavr)
print('R_avr_initial:',r_avr)
print ('Number of elements:', len(O.bodies))
Vt = triax.boxVolume
print ('Vt:', Vt)
Vs = triax.particlesVolume
print('Vs:',Vs)
Vv=Vt-Vs
print('Vv:',Vv)
n = triax.porosity
print ('n:',triax.porosity)
e=n/(1-n)
print ('e:',e)

#######################################
################# END #################

plt.show()

O.pause()