yade-users team mailing list archive

[ehsan benabbas <question688421@xxxxxxxxxxxxxxxxxxxxx>]

New question #688421 on Yade:
https://answers.launchpad.net/yade/+question/688421

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
from yade import pack

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

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

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

from yade.params import table

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=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)
#sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

from yade import export
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,
 stressMask = 7,
 internalCompaction=True,
)

################# 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 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
if nRead==0: yade.qt.Controller(), yade.qt.View()

################# 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 ('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[1]	###%%%***
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
	setContactFriction(radians(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 #################

RRmax=max([b.shape.radius for b in O.bodies])
RRmin=min([b.shape.radius for b in O.bodies])
print('Maximum Radius:',RRmax)
print('Minimum Radius:',RRmin)
print ('Number of elements:', len(O.bodies))
print ('Box Volume engine:', triax.boxVolume)
Vt = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])  # 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. ***")
print ('step that starts the deviatoric loading ', O.iter/ORN)




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 ##################
Maximum Radius: 0.0005368387128864783
Minimum Radius: 0
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
step that starts the deviatoric loading  635.0
============ APPLYING DEVIATORIC LOADING ============

Thank you,
Ehsan

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.