yade-users team mailing list archive

[Luc Scholtès <question222250@xxxxxxxxxxxxxxxxxxxxx>]

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

Luc Scholtès posted a new comment:
Hi Anna,

I could not fin an email adress where to send you my script and some
explanations, so here is my code below. I guess you can copy-paste it
directly as a python script. Be careful, some things are maybe
deprecated for last versions of YADE. TranslationEngine is used but, as
told by Bruno, it is not mandatory and could/should be replaced by a
condition on the boxes velocity. I am not exactly happy with the way the
simulation is servo-controlled as is, but it works if the loading
parameters are well adjusted (this is where I should take some time to
generalize the simulation like it is done for the triaxial test or
"Jerome's engine"). Be sure to adapt it to the interaction law you are
using in the definition of the "InteractionLoop" (I use here a modified
version of one of the currently available law).

Briefly, here is the process:
- packing is imported from a predefined text file. 
- a first box is appended at the top of the sample. 
- normal loading is performed until the predefined value is reached.
- once equilibrium is reached, top and bottom half boxes of a "Casagrande box" are created
- the bottom box is then displaced at constant velocity

Do not hesitate if you have any questions.

Cheers

  Luc

# -*- coding: utf-8 -*-
O=Omega()

from yade import ymport, utils, plot

################## parameters definition 
PACKING='X1Y1Z1_5k'
DAMPING=0.5
dtCOEFF=0.5
normalSTRESS=5e5
normalVEL=normalSTRESS/1e8 # 0.001 for 100kPa // optimized for normalVEL=normalSTRESS/1e8?
shearVEL=1*normalVEL # try different values to ensure quasi-static conditions
intR=1.263 #1.263 for X1Y1Z1_5k
DENS=3000
YOUNG=50e9
FRICT=18
ALPHA=0.3
TENS=45e5
COH=45e6
iterMAX=400000
Output='shearTest_'+PACKING+'_K10_D3000_Y50e9A03F18T45e5C45e6_500kPa_shearVel1'

###################### Import of the sphere assembly
### material definition
def sphereMat(): return CFpmMat(type=1,young=YOUNG,frictionAngle=radians(FRICT),density=DENS)
def wallMat(): return CFpmMat(type=0,young=YOUNG,frictionAngle=radians(0),density=DENS)

## copy spheres from the packing into the scene
O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

## preprocessing to get dimensions
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf
## initial surface
S0=X*Z

## spheres factory
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   o.shape.color=(0.7,0.5,0.3)
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/numSpheres

###################### creation of shear box
thickness=Y/100
oversizeFactor=1.3

### loading platens
O.bodies.append(utils.box(center=(xinf+X/2.,yinf-thickness+Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness,oversizeFactor*Z/2),material=wallMat,fixed=True)) 
bottomPlate=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2.,ysup+thickness-Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness,oversizeFactor*Z/2),material=wallMat,fixed=True))
topPlate=O.bodies[-1]

###################### Engines
O.engines=[
  
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='aabb'),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='ss2d3dg'),Ig2_Box_Sphere_ScGeom()],
		[Ip2_CFpmMat_CFpmMat_CFpmPhys(cohesiveTresholdIteration=1,alpha=ALPHA,tensileStrength=TENS,cohesion=COH,label='interactionPhys')],
		[Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM(recordCracks=recCRACKS,label='interactionLaw')]
	),

	GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),timestepSafetyCoefficient=dtCOEFF),
	TranslationEngine(ids=[topPlate.id],translationAxis=(0.,-1.,0.),velocity=0.,label='yTranslation'),
	PyRunner(iterPeriod=1,initRun=True,command='servoController()',label='servo'),
	NewtonIntegrator(damping=DAMPING,gravity=(0,0,0),label='damper'),
	PyRunner(iterPeriod=iterMAX/1000,initRun=True,command='dataCollector()'),
]

###################### Engines definition ( servoController() and dataCollector() )
shearing=False
sigmaN=0
tau=0
Fs1=0
Fs2=0
Xdispl=0
px0=0
Ydispl=0
py0=topPlate.state.pos[1]
prevTranslation=0
n=0

def servoController():
    global px0, py0, sigmaN, n, Fn1, Fn2, shearing, butee, piston1, piston2
    Fn1=abs(O.forces.f(topPlate.id)[1])
    Fn2=abs(O.forces.f(bottomPlate.id)[1])
    sigmaN=Fn1/(S0-2*Xdispl*Z) # necessary?
    #print yTranslation.velocity, sigmaN
    if shearing==False:
      if yTranslation.velocity<normalVEL:
	yTranslation.velocity+=normalVEL/1000
      if sigmaN>(0.9*normalSTRESS):
	yTranslation.velocity=normalVEL*((normalSTRESS-sigmaN)/normalSTRESS)
    if shearing==False and abs((normalSTRESS-sigmaN)/normalSTRESS)<0.001 and unbalancedForce()<0.01:
      yTranslation.velocity=0
      n+=1
      if n>1000 and abs((sigmaN-normalSTRESS)/normalSTRESS)<0.001:
	print 'stress on joint plane = ', utils.forcesOnPlane((X/2,Y/2,Z/2),(0,1,0))/S0
	### add top box
	O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,3*(ysup/4)+2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
	butee=O.bodies[-1]
	O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,3*(ysup/4)+2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
	O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))
	O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))                                                                                                                                                            
	### add bottom box
	O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,(ysup/4)-2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
	piston1=O.bodies[-1]
	O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,(ysup/4)-2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
	piston2=O.bodies[-1]
	O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))
	O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))
	### add engine and initialisation
	O.engines=O.engines[:6]+ [TranslationEngine(ids=[piston1.id,piston2.id],translationAxis=(-1.,0.,0.),velocity=0.)]+O.engines[6:]
	px0=piston1.state.pos[0]
	py0=topPlate.state.pos[1]
	shearing=True
	print 'shearing! || iteration=', O.iter
    if shearing==True:
      #yTranslation.velocity=0. # if you want a constant normal diaplacement control
      yTranslation.velocity=50*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
      if O.engines[6].velocity<shearVEL:
	O.engines[6].velocity+=(shearVEL/1000)

def dataCollector():
    global Xdispl, Ydispl, tau, Fs1, Fs2
    if shearing:
      Fs1=abs(O.forces.f(butee.id)[0])
      Fs2=abs(O.forces.f(piston1.id)[0])
      Xdispl=abs(piston1.state.pos[0]-px0)
    tau=Fs1/(S0-2*Xdispl*Z)
    Ydispl=abs(topPlate.state.pos[1]-py0)
    numberOfCohesiveBonds=0
    for i in O.interactions:
      if i.phys.isCohesive:
	numberOfCohesiveBonds+=1
    yade.plot.addData({'t':O.time,'i':O.iter,'Xdispl':Xdispl,'sigmaN':sigmaN,'Fn1':Fn1,'Fn2':Fn2,'tau':tau,'Fs1':Fs1,'Fs2':Fs2,'Ydispl':Ydispl,'unbF':utils.unbalancedForce(),'Ek':utils.kineticEnergy(),'bonds':numberOfCohesiveBonds})
    plot.saveGnuplot(Output)
    
# defines window plot
plot.plots={'i':('sigmaN','tau')}

################# graphical intervace
from yade import qt
qt.Controller()
qt.View()
v=qt.Renderer()
v.dispScale=(1,1,1) # displacements are scaled for visualization purpose

################# to manage interaction detection factor during the first timestep
O.step();
################# initializes the interaction detection factor to its default value (new contacts will be created between strictly contacting particles only)
ss2d3dg.interactionDetectionFactor=-1.
aabb.aabbEnlargeFactor=-1.

################# simulation starts here
O.run(iterMAX)

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.