New question #694100 on Yade:
https://answers.launchpad.net/yade/+question/694100
Goodmoring,
using this script [1] I've made some simulation with different value of SN
parameter and PI (25000, 50000 an 75000 for SN and respectively 2500, 5000 and
7500 for PI). By comparing those result with real shear test performed on a
Bristol Sand ( confining stress was set on 25, 50 and 75 kPa) I've observed
that the simulated stress-strain graph is exactly the half of the stress-strain
curve at each confining stress. I've tried to double SN and PI parameters on
the simulation, and the results o a stress-strain curve match the real
stress-strain curve at each confining stress. How is that possible? Is There
something I don't understand that makes this possible, like a "trick" on the
script?
Thanks for your help.
[1]
###
from yade import pack,plot,export
import math
sp=pack.SpherePack()
O.periodic=True
# dimensions of sample (fixed by particle size such as L/D~15)
RADIUS=0.001
length=0.06
height=length/3
width=length
thickness=RADIUS
# friction angles
compFRIC=15. # during compaction (controls porosity)
FRIC=40. # during shear
# boundary conditions
PI=5e3 # sample preparation: pressure applied for the isotropic compaction
SN=50e3 # normal stress
RATE=0.1 # shear velocity (top plate)
# simulation control
DAMPSHEAR=0
ITER=8e4
VTK=20
OUT='TEST'
####
O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=1600,young=230e6,poisson=0.3,frictionAngle=radians(compFRIC),normalCohesion=0,shearCohesion=0,label='sphereMat'))
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])
sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.001,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for
s in sp])
effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*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())
,PeriTriaxController(dynCell=True,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
,NewtonIntegrator(damping=0.3,label='newton')
,PyRunner(command='dataRecorder()',iterPeriod=100,label='recData',dead=False)
,VTKRecorder(fileName=OUT+'.',iterPeriod=1,skipNondynamic=1,recorders=['spheres','colors','velocity','bstresses','stress','boxes','intr','boxes'],label='saveSolid',dead=True)
]
def dataRecorder():
h=vol=vol_s=nb_s=0.
h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
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[0]!=1:
nb_s+=1
vol_s += 4.*pi/3.*(o.shape.radius)**3
n = 1-vol_s/vol
nbFrictCont=0.
for i in O.interactions:
if i.isReal and i.phys.cohesionBroken:
nbFrictCont+=1
plot.addData(
iter=O.iter
,contactStress01=((contactStress[0,1])/1000),contactStress11=(contactStress[1,1])
,xW=(O.bodies[0].state.pos[0]*1000)
,height=h
,volume=vol
,porosity=n
)
def triaxDone():
global phase
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
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+=o.shape.radius
if o.shape.radius>Rmax: Rmax=o.shape.radius
if o.shape.radius<Rmin: Rmin=o.shape.radius
vol_s += 4.*pi/3.*(o.shape.radius)**3
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)
triax.dead=True
O.pause()
#### Initialization
print 'SAMPLE PREPARATION!'
recData.dead=True # uncomment to record what is happening during stress
initialization
O.run(600000,1)
saveSolid.dead=False
#saveSolid.fileName=OUT+'_isoConfined.'
O.step()
saveSolid.dead=True
O.save(OUT+'_initialState.yade')
print 'Normal stress (platen) =
',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) =
',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]
#### Applying normal stress
print 'NORMAL LOADING! iter=',O.iter
recData.dead=True
stage=0
stiff=fnPlaten=currentSN=0.
def servo():
global stage,stiff,fnPlaten,currentSN
if stage==0:
currentSN=O.forces.f(0)[1]/(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[0].state.vel[1]=boundaryVel
if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
stage+=1
fnPlaten=O.forces.f(0)[1]
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[0].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)[1]-fnDesired)/stiff/O.dt),O.forces.f(0)[1]-fnDesired)
boundaryVel=copysign(min(RATE,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
O.bodies[0].state.vel[1]=boundaryVel
O.engines =
O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]
O.run(35000,1)
print 'STABILIZING! iter=',O.iter
O.run(10000,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[0]<xmin: xmin=o.state.pos[0]
if (o.state.pos[0]>=i*dxi) and (o.state.pos[0]<((i+0.5)*dxi)):
o.shape.color[1]=1
for o in O.bodies:
if (o.state.pos[0]>(xmin+0.5*O.cell.hSize[0,0]-RADIUS)) and
(o.state.pos[0]<(xmin+0.5*O.cell.hSize[0,0]+3*RADIUS)):
o.shape.color[2]=0
saveSolid.dead=False
#saveSolid.fileName=OUT+'_normalLoaded.'
O.step()
saveSolid.dead=True
O.save(OUT+'_normallyLoaded.yade')
if recData.dead==False: plot.saveDataTxt(OUT)
print 'Normal stress (platen) =
',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) =
',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*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[0]=1
gluedSpheres += [O.bodies[i.id2]]
if isinstance(O.bodies[i.id2].shape,Box):
O.bodies[i.id1].shape.color[0]=1
gluedSpheres += [O.bodies[i.id1]]
for i in O.interactions:
if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and
O.bodies[i.id2].shape.color[0]==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
saveSolid.dead=False
O.step()
saveSolid.dead=True
O.save(OUT+'_shearInit.yade')
recData.dead=False
newton.damping=DAMPSHEAR
saveSolid.dead=False
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[0].state.vel[0]=shearVel
## top and bottom walls move
#O.bodies[0].state.vel[0]=0.5*shearVel
#O.bodies[1].state.vel[0]=-0.5*shearVel
if ( O.iter >= (int(iterShear+ITER)) ):
print 'iter=',O.iter,' -> FINISHED!'
plot.saveDataTxt(OUT)
O.save(OUT+'_sheared.yade')
sys.exit(0)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
