yade-users team mailing list archive

[Luc Scholtès <question672563@xxxxxxxxxxxxxxxxxxxxx>]

Question #672563 on Yade changed:
https://answers.launchpad.net/yade/+question/672563

Luc Scholtès proposed the following answer:
In my mind, the biperiodic configuration where the sample is sandwiched
in between 2 moving plates corresponds to a direct shear test... It
does, right? Also, I guess that, instead of a translation, you could try
to apply a rotation to the top plate with the centre of rotation located
either the centre of the sample or offset from it so as to mimic your
experiment...

ANW, here is the code for simulating the "sandwiched sample". It is a
bit long but it works and produces outputs that should be useful.

Cheers

###

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.05
length=15*(2*RADIUS)
height=length 
width=length
thickness=RADIUS

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

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=2500,young=1e8,poisson=0.5,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.2,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,mass=10,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=10,label='recData',dead=True)
	,VTKRecorder(fileName=OUT+'.',iterPeriod=1,skipNondynamic=1,recorders=['spheres','colors','velocity','bstresses'],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
		,stress_upWall0=abs(O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall1=abs(O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall2=abs(O.forces.f(0)[2]/(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[0].state.pos[0]
		,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[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=False # uncomment to record what is happening during stress initialization 
O.run(1000000,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

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(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[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.