Looking for Direct Shear test script

Asked by Amedeo Galletti on 2020-10-05

I was looking for a direct shear test script on https://yade-dem.org/doc/tutorial.html , I only found https://yade-dem.org/doc/tutorial-examples.html#periodic-simple-shear . Do you know where I can find a script that simulate a direct shear test?

Thanks

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-10-05
Last reply:
2020-10-12
Jan Stránský (honzik) said : #1

Hello,

according to the amount of answers, the question provided is not motivating enough to answer.
Please provide more information, for example provide (not only!):
- material
   - real - soil, concrete, ...
   - how you want to simulate it
- 2D / 3D
- loading
   - lateral stress?
   - ...
- what results you are interested in
   - forces
   - failure mode
   - ...
- do you have some ideas / have you tried something already?
   - e.g. boundary conditions - using walls or directly prescribed motion of boundary spheres or ...
   - ...
- what have you searched (googling 'yade "direct shear"' or 'dem "direct shear"' gives plenty of results) and why it did not help?
- ...
- ...

cheers
Jan

Luc Scholtès (luc) said : #2

Hi,

The following lines may help you build a direct shear test simulation (the script needs other scripts to work properly, i.e., one to generate the sample and one to do the identification of the contact along the joint) . Let us know if that's what you are after.

###

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

from yade import ymport, utils, plot

################## parameters definition
PACKING='global_spheres.xyz' # imported from a text file

intR=1.4
DENS=3000
YOUNG=12e9
FRICT=18
POISSON=0.15
TENS=11e6
COH=11e7
SMOOTH=True
jointNSTIFF=12e9
jointSSTIFF=12e8
jointFRICT=30
jointTENS=0.
jointCOH=0.
jointDIL=0.

DAMPING=0.5
dtCOEFF=0.5
iterMAX=100000
saveVTK=10
recCRACKS=True
Output=PACKING+'_1MPa'

normalSTRESS=1.e6
normalVEL=normalSTRESS/1e8 # optimized for normalVEL=normalSTRESS/1e8?
shearVEL=2*normalVEL # optimized for shearVEL=normalVEL?

###################### Import of the sphere assembly
### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=POISSON,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH,jointNormalStiffness=jointNSTIFF,jointShearStiffness=jointSSTIFF,jointTensileStrength=jointTENS,jointCohesion=jointCOH,jointFrictionAngle=radians(jointFRICT),jointDilationAngle=radians(jointDIL))
def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,poisson=POISSON,frictionAngle=radians(0))

## copy spheres from the packing into the scene
O.bodies.append(ymport.text(PACKING,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

#### create the fracture directly here
import gts
ptA = gts.Vertex( xinf-X/4., yinf-Y/4., zinf+Z/2.)
ptB = gts.Vertex( xinf-X/4., ysup+Y/4., zinf+Z/2.)
ptC = gts.Vertex( xsup+X/4., yinf-Y/4., zinf+Z/2.)
ptD = gts.Vertex( xsup+X/4., ysup+Y/4., zinf+Z/2.)
e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptC)
e3 = gts.Edge(ptC,ptB)
f1 = gts.Face(e1,e2,e3)
e4 = gts.Edge(ptC,ptD)
e5 = gts.Edge(ptD,ptB)
f2 = gts.Face(e4,e5,e3)
s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)

execfile('identifBis.py') #you need the script in the same folder (this one is available in the JCFPM examples)

###################### creation of shear box
thickness=Z/100.
oversizeFactor=1.3

### loading platens
O.bodies.append(utils.box(center=(xinf+X/2.,yinf+Y/2.,zinf-thickness+Rmean/5.),extents=(oversizeFactor*X/2,oversizeFactor*Y/2,thickness),material=wallMat,fixed=True))
bottomPlate=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2.,yinf+Y/2.,zsup+thickness-Rmean/5.),extents=(oversizeFactor*X/2,oversizeFactor*Y/2,thickness),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_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=SMOOTH,recordCracks=True,label='interactionLaw')]
 ),

 GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),timestepSafetyCoefficient=dtCOEFF),
 TranslationEngine(ids=[topPlate.id],translationAxis=(0.,0.,-1.),velocity=0.,label='zTranslation'),
 PyRunner(iterPeriod=1,initRun=True,command='servoController()',label='servo'),
 NewtonIntegrator(damping=DAMPING,gravity=(0,0,0),label='damper'),
 PyRunner(iterPeriod=100,initRun=True,command='dataCollector()'),
 VTKRecorder(iterPeriod=int(iterMAX/saveVTK),initRun=True,fileName=Output+'-',recorders=['spheres','stress','velocity','jcfpm','cracks','boxes'])

]

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

def servoController():
    global px0, pz0, sigmaN, n, Fn1, Fn2, shearing, butee, piston1, piston2
    Fn1=abs(O.forces.f(topPlate.id)[2])
    Fn2=abs(O.forces.f(bottomPlate.id)[2])
    sigmaN=Fn1/(S0-2*Xdispl*Y) # necessary?
    #print yTranslation.velocity, sigmaN
    if shearing==False:
      if zTranslation.velocity<normalVEL:
 zTranslation.velocity+=normalVEL/100.
      if sigmaN>(0.95*normalSTRESS):
 zTranslation.velocity=10*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
    if shearing==False and abs((normalSTRESS-sigmaN)/normalSTRESS)<0.001 and unbalancedForce()<0.005:
      zTranslation.velocity=0
      print 'stress on joint plane = ', utils.forcesOnPlane((X/2,Y/2,Z/2),(0,0,1))/S0
      ### top box
      O.bodies.append(utils.box(center=(xinf-thickness+Rmean/5.,(yinf+0.5*Y),(zinf+0.75*Z)+3*Rmean)
    ,extents=(thickness,oversizeFactor*Y/2.,Z/4.),material=wallMat,fixed=True))
      butee=O.bodies[-1]
      O.bodies.append(utils.box(center=(xsup+thickness-Rmean/5.,(yinf+0.5*Y),(zinf+0.75*Z)+3*Rmean),
    extents=(thickness,oversizeFactor*Z/2.,Z/4.),material=wallMat,fixed=True))
      O.bodies.append(utils.box(center=((xinf+0.5*X),yinf-thickness+Rmean/5.,(zinf+0.75*Z)+3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4.),material=wallMat,fixed=True,wire=True))
      O.bodies.append(utils.box(center=((xinf+0.5*X),ysup+thickness-Rmean/5.,(zinf+0.75*Z)+3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4.),material=wallMat,fixed=True,wire=True))
      ### bottom box
      O.bodies.append(utils.box(center=(xsup+thickness-Rmean/5.,(yinf+0.5*Y),(zinf+0.25*Z)-3*Rmean),
    extents=(thickness,oversizeFactor*Y/2.,Z/4.),material=wallMat,fixed=True))
      piston1=O.bodies[-1]
      O.bodies.append(utils.box(center=(xinf-thickness+Rmean/5.,(yinf+0.5*Y),(zinf+0.25*Z)-3*Rmean),
    extents=(thickness,oversizeFactor*Y/2.,Z/4.),material=wallMat,fixed=True))
      piston2=O.bodies[-1]
      O.bodies.append(utils.box(center=((xinf+0.5*X),yinf-thickness+Rmean/5.,(zinf+0.25*Z)-3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4),material=wallMat,fixed=True,wire=True))
      O.bodies.append(utils.box(center=((xinf+0.5*X),ysup+thickness-Rmean/5.,(zinf+0.25*Z)-3*Rmean),
    extents=(oversizeFactor*X/2.,thickness,Z/4.),material=wallMat,fixed=True,wire=True))
      ### initialisation
      O.engines=O.engines[:5]+ [TranslationEngine(ids=[piston1.id,piston2.id],translationAxis=(-1.,0.,0.),velocity=0.)]+O.engines[5:]
      px0=piston1.state.pos[0]
      pz0=topPlate.state.pos[2]
      shearing=True
      print 'shearing now! || iteration=', O.iter
    if shearing==True:
      zTranslation.velocity=0.
      ## Removes collider (lattice like) -> not really helpful in this particular case
      #O.engines=O.engines[:1]+O.engines[2:]
      zTranslation.velocity=50*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
      if O.engines[5].velocity<shearVEL:
 O.engines[5].velocity+=(shearVEL/1000)

def dataCollector():
    global Xdispl, Zdispl, 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*Y) # necessary?
    Zdispl=abs(topPlate.state.pos[2]-pz0)
    yade.plot.addData({'t':O.time,'i':O.iter,'Xdispl':Xdispl,'sigmaN':sigmaN,'tau':tau,'Zdispl':Zdispl,'unbF':utils.unbalancedForce(),'tc':interactionLaw.nbTensCracks,'sc':interactionLaw.nbShearCracks})
    plot.saveDataTxt(Output)
    #plot.saveGnuplot(Output)

################# 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(int(iterMAX))

Amedeo Galletti (amedeog) said : #3

Hello,
Sorry for the lack of infos. I’m new to the yade community.

Luc, thanks for your support. I’ve used a script that you posted on this forum:

###

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)

Is it suitable for a model of a direct shear test?
What about units of measurement? What’s the unit of RADIUS for example.

Thanks, and sorry for my poor English.

Amedeo

Jan Stránský (honzik) said : #4

> What about units of measurement? What’s the unit of RADIUS for example.

Yade has no units, just numbers [1,2,3,4].
If you use basic SI units for all inputs, you get all output in basic SI units, too.
Using only basic SI units is a good practice.

cheers
Jan

[1] https://answers.launchpad.net/yade/+question/691104
[2] https://answers.launchpad.net/yade/+question/690022
[3] https://answers.launchpad.net/yade/+question/687871
[4] https://answers.launchpad.net/yade/+question/672270

Can you help with this problem?

Provide an answer of your own, or ask Amedeo Galletti for more information if necessary.

To post a message you must log in.