Modeling a Direct Shear Test with Simple Shear

Hi,
You may already figured out a solution for your problem but here is nonetheless a suggestion. You could create a simulation where a (bi)periodic packing is sandwiched between two horizontal plates. The packing would thus have periodic boundaries along the 2 horizontal directions while the upper and bottom (infinite) plates will allow you to apply a given normal stress on it and to shear it. Once the desired normal stress achieved, move the plates to produce shear. All this can be achieve, e.g., by fixing the bottom plate and by servo-controlling the upper plate (vertical displacement to control the normal stress and horizontal displacement to shear the sample). Depending on your need, you may increase friction between the plates and the sample by fixing the contacting particles to the plates.
I can provide a script if you think that's appropriate.
Luc

Luc,

Thank you for responding. A mixed boundary condition environment is a great
idea, and it may be the route that I end up taking. If it's not too much
trouble, could you send the example script that you offered? This would
help out tremendously with possible troubleshooting.

The solution you proposed is based on simple shear. Do you think it's
appropriate to model direct/rotational shear by simplifying it to a simple
shear scenario? Let me know what you think.

Thanks,
Austin

On Thu, Aug 30, 2018 at 4:02 AM, Luc Scholtès <
<email address hidden>> wrote:

> Your question #672563 on Yade changed:
> https://answers.launchpad.net/yade/+question/672563
>
> Status: Open => Answered
>
> Luc Scholtès proposed the following answer:
> Hi,
> You may already figured out a solution for your problem but here is
> nonetheless a suggestion. You could create a simulation where a
> (bi)periodic packing is sandwiched between two horizontal plates. The
> packing would thus have periodic boundaries along the 2 horizontal
> directions while the upper and bottom (infinite) plates will allow you to
> apply a given normal stress on it and to shear it. Once the desired normal
> stress achieved, move the plates to produce shear. All this can be achieve,
> e.g., by fixing the bottom plate and by servo-controlling the upper plate
> (vertical displacement to control the normal stress and horizontal
> displacement to shear the sample). Depending on your need, you may increase
> friction between the plates and the sample by fixing the contacting
> particles to the plates.
> I can provide a script if you think that's appropriate.
> Luc
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/672563/+confirm?answer_id=0
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/672563
>
> You received this question notification because you asked the question.
>

--
*Austin Sutton*
Graduate Research Assistant and Graduate Teaching Assistant
342 Toomey Hall, Rolla, MO 65409
PhD Student - Mechanical Engineering Department
Missouri University of Science and Technology

>Do you think it's appropriate to model direct/rotational shear by simplifying it to a simple
shear scenario?

Hi,
They are fundamentally different I'm afraid. Even direct and rotational are not the same.
All three types of boundary conditions can be (and have been) simulated with DEM.

There is no general and convincing argument to say they are all equivalent (instead, I've always heard that they were different...).
I guess finding justifications (if any) will depend on your particular context.
Sorry for late reply
Bruno

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)

>the biperiodic configuration where the sample is sandwiched in between 2 moving plates corresponds to a direct shear test...

Direct shear refers to the (inhomogeneous) shear box for soils, or a sheared rock joint.
The yade's preprocessor "SimpleShear" is badly named, btw, since it simulates a direct shear.
Rotational shear is not homogeneous, although it does not impose the position of the shear band as firmly as direct shear.

Conversely, the periodic BCs define simple shear (homogeneous). Only very specific devices (e.g. 1γ2ε in Grenoble) can do such thing in practice.

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

Nope. Such a rotational field is not periodic, the simulation would just go weird. It is not possible to use periodic BCs for rotational shear (unless someone implements periodicity for annular geometries, which is not straightforward).

Bruno

*Luc,* I appreciate your effort in putting together that script. It is
long, but it I'm grateful since it provides a lot of details that I was
going to need to figure out anyways. Thanks!

--> 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?
After reading a bit more into the topic, it seems that it is simple shear.
However, the script you made is still useful since I've read that direct
and simple shear tests have been shown to provide close results in terms of
steady-state stress for some materials. In other words, they may differ in
terms of peak stress and strain, but can give the same yield locus curves.

*Bruno,* thank you for your input.

--> They are fundamentally different I'm afraid. Even direct and rotational
are not the same.
This is exactly what I was thinking, but wasn't entirely sure. Thank you
for the clarification.

--> The yade's preprocessor "SimpleShear" is badly named, btw, since it
simulates a direct shear.
I had no idea. I dismissed this route because of its name...thank you for
bringing this to my attention.

*To all,* thanks for your input. I believe I have a couple of different
routes to try. One last question: have either of you heard of reducing the
modulus of the powder material in a simulated shear cell to reduce
computational time? I'm just wondering since my material is small stainless
steel particles: high modulus with a minimum particle size of 1 um == very
small timestep. Since we care about stresses, I'm wondering if this would
be frowned upon. If I were to reduce the modulus, I believe I would just
need to ensure that the volume ratio of the packing when initially
consolidated is the same as that observed experimentally.

On Fri, Aug 31, 2018 at 10:43 AM, Bruno Chareyre <
<email address hidden>> wrote:

> Your question #672563 on Yade changed:
> https://answers.launchpad.net/yade/+question/672563
>
> Bruno Chareyre proposed the following answer:
> >the biperiodic configuration where the sample is sandwiched in between
> 2 moving plates corresponds to a direct shear test...
>
> Direct shear refers to the (inhomogeneous) shear box for soils, or a
> sheared rock joint.
> The yade's preprocessor "SimpleShear" is badly named, btw, since it
> simulates a direct shear.
> Rotational shear is not homogeneous, although it does not impose the
> position of the shear band as firmly as direct shear.
>
> Conversely, the periodic BCs define simple shear (homogeneous). Only
> very specific devices (e.g. 1γ2ε in Grenoble) can do such thing in
> practice.
>
> > 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...
>
> Nope. Such a rotational field is not periodic, the simulation would just
> go weird. It is not possible to use periodic BCs for rotational shear
> (unless someone implements periodicity for annular geometries, which is
> not straightforward).
>
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/672563/+confirm?answer_id=4
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/672563
>
> You received this question notification because you asked the question.
>

--
*Austin Sutton*
Graduate Research Assistant and Graduate Teaching Assistant
342 Toomey Hall, Rolla, MO 65409
PhD Student - Mechanical Engineering Department
Missouri University of Science and Technology

@Bruno (and everyone interested in this topic) :

The box defined by the preprocessor SimpleShear can be sheared as depicted Fig. 2 (right) of [*] (provided the user adds KinemCNDEngine [**] to O.engines)

To me, this looks exactly the same than a simple shear box like the 1γ2ε, and quite different than a direct shear box.

[*] https://www.sciencedirect.com/science/article/pii/S1365160910001784
[**] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.KinemCNDEngine

>The box defined by the preprocessor SimpleShear can be sheared as depicted Fig. 2 (right) of [*] (provided the user adds KinemCNDEngine [**] to O.engines)

Hi Jérôme,
Then it is my mistake, sorry.
B

p.s. I don't understand "can be ... provided the user adds". A preprocessor is supposed to be ready to run, not the case here?
Also note that it is broken (FATAL /data/trunk/core/ThreadRunner.cpp:30 run: Exception occured: output stream error).

The preprocessor should run, and "carries out an oedometric compression, until a value of normal stress equal to 2 MPa (and a stable mechanical state)".

The fact is there is no proper shear simulated with this preprocessor (like if TriaxialTest would just do the isotropic compression stage..), in this sense it's true it's badly named ;-)
The doc [*] says it though.

As for the crash, thanks for reporting. However, here (with yadedaily 2018.02b-1727253714~xenial, and trunk 2018-08-06.git-7253714) it works (?)

[*] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.SimpleShear

This question was expired because it remained in the 'Open' state without activity for the last 15 days.