# Relative Density

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2020-10-20: #1

Hello,

> How can I ...

- Why do you want to control it?
- At what circumstances (initial packing, during simulation)?
- ...

> control

what does "control" mean? check (get)? prescribe?

> in a script

do you have any specific script, or is it just a general question?

> Is relDensity good?

depends on "control", see above

cheers
Jan

 Revision history for this message Amedeo Galletti (amedeog) said on 2020-10-21: #2

###

import math

sp=pack.SpherePack()
O.periodic=True

# dimensions of sample (fixed by particle size such as L/D~15)
height=length
width=length

# friction angles
compFRIC=10. # during compaction (controls porosity)
FRIC=30. # during shear

# boundary conditions
PI=1.e5 # sample preparation: pressure applied for the isotropic compaction
SN=5.e6 # normal stress
RATE=0.1 # shear velocity (top plate)

# simulation control
DAMPSHEAR=0.
ITER=2e5
VTK=20
OUT='TEST'

####

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

upBox = utils.box(center=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])

O.bodies.append([utils.sphere(s,s,color=(0,0,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies.state.pos-O.bodies.state.pos)*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

#print 'volRatio=',volRatio

O.engines=[
ForceResetter()
,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
,InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
)
,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
,NewtonIntegrator(damping=0.3,label='newton')
]

def dataRecorder():
h=vol=vol_s=nb_s=0.
h=O.bodies.state.pos-O.bodies.state.pos
vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
contactStress=getStress(vol)
for o in O.bodies:
if isinstance(o.shape,Sphere) and o.shape.color!=1:
nb_s+=1
n = 1-vol_s/vol
nbFrictCont=0.
for i in O.interactions:
if i.isReal and i.phys.cohesionBroken:
nbFrictCont+=1
iter=O.iter
,stress_upWall0=abs(O.forces.f(0)/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall1=abs(O.forces.f(0)/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall2=abs(O.forces.f(0)/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
,contactStress00=(contactStress[0,0]),contactStress01=(contactStress[0,1]),contactStress02=(contactStress[0,2]),contactStress10=(contactStress[1,0]),contactStress11=(contactStress[1,1]),contactStress12=(contactStress[1,2]),contactStress20=(contactStress[2,0]),contactStress21=(contactStress[2,1]),contactStress22=(contactStress[2,2])
,xW=O.bodies.state.pos
,height=h
,volume=vol
,porosity=n
,k=2.0*nbFrictCont/nb_s
,unbF=unbalancedForce()
)

def triaxDone():
global phase
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies.state.pos-O.bodies.state.pos)*O.cell.hSize[0,0]*O.cell.hSize[2,2])
h=O.bodies.state.pos-O.bodies.state.pos
vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
contactStress=getStress(vol)
vol_s=Rmean=Rmax=nbSph=0
Rmin=1e6
for o in O.bodies:
if isinstance(o.shape,Sphere):
nbSph+=1
Rmean=Rmean/nbSph
n = 1-vol_s/vol
print 'DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n
print 'Changing contact properties now'
tt=TriaxialCompressionEngine()
tt.setContactProperties(FRIC)
O.pause()

#### Initialization
print 'SAMPLE PREPARATION!'

#recData.dead=False # uncomment to record what is happening during stress initialization
O.run(1000000,1)
#saveSolid.fileName=OUT+'_isoConfined.'
O.step()

print 'Normal stress (platen) = ',O.forces.f(0)/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies.state.pos-O.bodies.state.pos)*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### Applying normal stress

stage=0
stiff=fnPlaten=currentSN=0.
def servo():
global stage,stiff,fnPlaten,currentSN
if stage==0:
currentSN=O.forces.f(0)/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
unbF=unbalancedForce()
#print 'SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF
boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
O.bodies.state.vel=boundaryVel
if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
stage+=1
fnPlaten=O.forces.f(0)
print 'Normal stress =',currentSN,' | unbF=',unbF
## the following computes the stiffness of the plate (used for stress control of the top plate)
for i in O.interactions.withBody(O.bodies.id):
stiff+=i.phys.kn
print 'DONE! iter=',O.iter
O.pause()
if stage==1:
fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
#boundaryVel=copysign(abs(0.5*(O.forces.f(0)-fnDesired)/stiff/O.dt),O.forces.f(0)-fnDesired)
boundaryVel=copysign(min(RATE,abs(0.333*(O.forces.f(0)-fnDesired)/stiff/O.dt)),O.forces.f(0)-fnDesired)
O.bodies.state.vel=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]

O.run(1000000,1)
print 'STABILIZING! iter=',O.iter
O.run(1000,1)
# coloring particles to see vertical stripes in material
dxi=4*(2*RADIUS) # can be a function of cell length dxi=O.cell.hSize[0,0]/5.
n=int(length/dxi)
xmin=1e6
for i in range(0,n):
for o in O.bodies:
if o.id>1:
if o.state.pos<xmin: xmin=o.state.pos
if (o.state.pos>=i*dxi) and (o.state.pos<((i+0.5)*dxi)):
o.shape.color=1
for o in O.bodies:
o.shape.color=0

O.step()

print 'Normal stress (platen) = ',O.forces.f(0)/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies.state.pos-O.bodies.state.pos)*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### preparing for shearing
Gl1_Sphere.stripes=1
print 'Gluing spheres to boundary walls'
gluedSpheres=[]

## gluing particles in contact with the walls
for i in O.interactions:
if i.isReal:
if isinstance(O.bodies[i.id1].shape,Box):
O.bodies[i.id2].shape.color=1
gluedSpheres += [O.bodies[i.id2]]
if isinstance(O.bodies[i.id2].shape,Box):
O.bodies[i.id1].shape.color=1
gluedSpheres += [O.bodies[i.id1]]

for i in O.interactions:
if i.isReal and ( O.bodies[i.id1].shape.color==1 and O.bodies[i.id2].shape.color==1 ):
O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
i.phys.initCohesion=True

print 'nb of glued spheres=',len(gluedSpheres)

#### shearing
print 'SHEARING! iter=',O.iter

O.step()

newton.damping=DAMPSHEAR
saveSolid.iterPeriod=int(ITER/VTK)
shearVel=0
iterShear=O.iter
while 1:
O.run(int(10),1)
if shearVel<RATE:
shearVel+=(RATE/100.)
# only top wall moves
O.bodies.state.vel=shearVel
## top and bottom walls move
#O.bodies.state.vel=0.5*shearVel
#O.bodies.state.vel=-0.5*shearVel
if ( O.iter >= (int(iterShear+ITER)) ):
print 'iter=',O.iter,' -> FINISHED!'
plot.saveDataTxt(OUT)
sys.exit(0)

Im refering to this script. I need to control the Relative Density of the sample.
For example for a Sand of Dr=40% .

 Revision history for this message Jan Stránský (honzik) said on 2020-10-21: #3

Hello,

> How can I ...

- Why do you want to control it?
- At what circumstances (initial packing, during simulation)?
- ...

> control

what does "control" mean? check (get)? prescribe?

> Is relDensity good?

depends on "control", see above

cheers
Jan

PS: is this message nearly the same as my first post? Yes, because you did not answer most of the questions. Posting 240 lines of code really does not help much. At least try to make it MWE, M=minimal .
One of the question is what "control" does mean. Your answer: "I need to control..."

 Revision history for this message Amedeo Galletti (amedeog) said on 2020-10-21: #4

Hi,
I want to create a pack of spheres with a certain relative density (for example a Sand with Dr=40%) . I want to do it at the beginning (when I'm creating the sample). My goal is to recreate a direct shear test of a sand, I think that Dr is an important parameter in order to recreate a sand.

I also found this on the Yade documentation , but I'm not able to fit it into the script that I posted earlier. I'm even not sure that is the right tool to achieve Relative Density of my sample.

Anyway thanks for your patience. I hope that I explained a bit better.

Regards,
Amedeo

 Revision history for this message Jan Stránský (honzik) said on 2020-10-21: #5

Thanks, much more understandable.

SpherePack.relDensity, as the docs suggests, is a "getter" - it returns computed value of relative density of the packing
So yes, it can be good for checking the value.

Since you are using PeriTriaxController, you can use it as strain-controlled (not stress controlled) compression to easily set target dimensions (knowing total volume of spheres) to get the desired relative density / porosity / void ratio ...

cheers
Jan