model a direct shear test

Hi Anna,

Just modify the example of the simple shear test ( in the Example folder of YADE). All required engines are already there.

The example of direct shear test can be included in the example part of YADE. You can write it totally in Python or generate one C++ engine and one simpler python script.

Viet

Hi Viet,
Thanks for your answer.
The problem is that the engine used in the simple shear test, if I have well understood, includes a rotation of the side walls and I am interested in an horizontal displacement. I am interested in the same kind of work as the one of Luc Scholtès describes in the paper Modelling progressive failure in fractured rock masses using a 3D discrete element method.
Anna

Hi,

I think, you do not need to use a C++ to solve this task.
You need to define to boxes (or set of facets) and use 2 TranslationEngines to
move them. You can also control their movement using engine labels and PyRunner.

==========================================
...
TranslationEngine(translationAxis=[1,0,0],velocity=2.0,ids=UpperBoxIDs,
label="Upper"),
TranslationEngine(translationAxis=[-1,0,0],velocity=2.0,ids=BottomBoxIDs,
label="Bottom"),
PyRunner(iterPeriod=50,command='ControlPlates()'),
...

def ControlPlates()
if (...):
  Upper.velocity = Upper.velocity - 1.0
  Bottom.velocity = Bottom.velocity - 1.0
==========================================

Cheers,

Anton

Hello Anna,

sorry for the late reply. I saw your email but could not find time to answer it.

Euh... Me again. I was looking for a way to attach files to my answer...

I have a python script that should help you. It performs a direct shear test as you describe with, first, the application of the normal load and, second, the shearing phase.

The problem is that, I am not very happy with this script as the control of the normal load during shearing is not "general". I wanted to improve it so as to do something more "general" like the triaxial engine for example but did not have the time. Currently, I have to adjust a variable that I have called servoCoeff depending on the parameters of the simulation (normal load, strength of the matrial, etc...).

Anyway, maybe we can discuss this by email so that I can send you the script (unless I can do it through launchpad?).

Cheers

Luc

Dear Anna,
Please allow me a question back. How did you have this idea of using translation engine? You are not alone. Actually it seems 99% of the users will jump on TranslationEngine as soon as they want to move an object, so we endlessly repeat that it is not needed. There is something to improve in the documentation or example scripts apparently, but I fail to identify what is misleading so many users.

Hi all,

@ Anton: I did it in a similar way but was not really sure.

@Luc: Of course, we can discuss this by email.

@Bruno: I don't know for me it was logical to use the TranslationEngine in order to translate the boxes. Could please tell me how you will do it without the translationEngine? Maybe just by changing directly the values of the velocity of the boxes?

cheers,
Anna

>Maybe just by changing directly the values of the velocity of the boxes?

Exactly. It is more direct (hence more predictable) and shorter to type.
How did you know a TranslationEngine exists in Yade? Reading the whole documentation?

Hello Anna,

I performed similar simulations, but with a simple shear box. I wrote
an engine to apply this constant normal stress path :
https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.KinemCNLEngine

There is an example of "my" simulations in
/examples/simple-shear/simpleShear.py

Our boxes are different (and I'm at the moment quite far from Yade),
but maybe this could help you..

Jerome

Hi guys,

Thank you for your help.

cheers,
Anna

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)

Hi Luc Scholtès,
Something is wrong when it go to the sentence O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
can you give me some hints how to build the PackingSpheres , I want to test it. thank you!

Hello,

Nothing is "wrong" with the code, it is simply looking for a packing whose name is defined by PACKING at the beginning of the script and finished by ".spheres". They are a lot of different ways to generate a packing (see: https://yade-dem.org/doc/user.html#sphere-packings). Usually, I first generate packings separately and I export tehm in a text file "name.spheres" so that I can use them for whatever I want later on. So you need to figure out how to generate your textfile.spheres using the method that better suit you using the utils function export.text. Another way, if you run a triaxial test using Yade graphical interface (ype F12 once yade terminal is open), a textfile.spheres is generated with predefined geometrical properties. You can work this way for parallellepipedic packings if you want.

Cheers

Luc