[Question #695561]: Make Radius Ratio for the sample and verity of particles radii

```New question #695561 on Yade:

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?

######################################################################################################
######### 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 pprint import pprint
O.timingEnabled=True

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

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

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

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
strainRate=-0.5  # %/sec     ### FIXME
target_strain=0.4 ### FIXME   %
sigmaIso=-5e5 # Pa=N/m2     ### FIXME
particleDensity=2600

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

print ('============ DEFINING FUNCTIONS ============')
def history():
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=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)

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

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/young,
thickness = thick,
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()]),
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=IPmacro,command='history()',label='macro_recorder'),
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
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 #################

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)

######################################################

triax.internalCompaction=False

# 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.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 ######################

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

--