How to make sure that the interaction between spheres and containers is frictionless

Hi,

Even though this would only add to the questions you have on the mailing list ;-), it would be better to make different Launchpad questions with these different questions :-)

1. and 3. and general answer : YADE bodies' shape is handled by Ig2 functors. Provided you can
- apply the same Material to your sphere and "container" bodies (possibly with different parameters values)
- and define a unique interaction geometrical description (like penetrationDepth which can be defined both for sphere-sphere, or sphere-facet interactions) to store in a unique IPhys

=> you do not need several (Ip2, Law2) functors, but just one pair (Ip2, Law2) to deal with sphere-sphere and sphere-container interactions.

Triaxial example script [1] for instance uses zero friction for the container bodies, and non-zero friction bodies for the sphere bodies. According to the logic implemented in Ip2_FrictMat_FrictMat_FrictPhys(), this leads to sphere-container interactions with a completely zero friction.

If you're using your own models / c++ classes, I would say it's up to you to imagine and code a suitable logic but it should be possible.

This being said, did I get correctly that you would like to have frictionless interactions at the boundaries, but still measure shear stress on the boundaries ? This is not possible.

2. Do you mean your shear box get open at one point ? If yes, you can not avoid problems. Would you leave sand flowing outside a direct shear box in the lab ? ;-)

Jérôme

[1] https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py

Hi Jerome,

Thanks so much. I checked the script you sent me.

1- I don't understand what is the meaning of giving same material parameters to both container and spheres. It should be logically wrong if you assign the same young and same Poisson to both when we know that the stiffness of steel is much higher than that of the soil. Also in my contact model Poisson is the ratio of the normal to the tangential contact stiffness (for my soil is 1.5 and I can't use the same value for steel-sphere contact)

2- I still am not able to use only one IP2 for both sphere-sphere and sphere-container since my Ip2 for sphere-sphere has some model specific parameters (I'll attach a part of my script for clarification) which are not applicable to sphere-container interaction.

3- I still need an absolute frictionless contact between spheres and containers. As you know in direct shear we really don't measure shear (frictional force between particles and some part of containers). In my case after reaching a target normal stress I move the lower box in horizontal direction (-X direction) and measure the normal force in the same direction on some parts of the upper box. This is the normal force on that plate (-X direction) induced by shearing the soil.

4- Yes, the shear box gets open at one point and this opening continues. If it is aimed to measure the post peak response it is necessary to continue shearing more and more up to some well published shear strain (let’s say 13%). As I move the lower box there will be two ever growing (I mean the area) open surfaces on the right and on the left. On the left, the bottom of some particles in the upper box are exposed to gravity (I have considered zero gravity but they still fall) and on the right the top of some particles in the lower box are exposed to free surface. In experimental direct shear this doesn't happen up to a certain horizontal displacement since shear box in the lab has thickness which avoids both issues. I have considered some plates to resemble the thickness of upper and lower shear boxes but after the shear box gets open the sphere come to contact with these plates and here is exactly where I need absolute frictionless contact otherwise the shear (friction) between particles and plates. Induce additional shear force which is not realistic since in the lab the contact between the thickness of shear box and spheres is completely frictionless. More importantly is the opening of lower shear box on the left. Since the particles have experienced some normal stress during consolidation they store some overlaps and when they find the opportunity to release this (exactly when the movement of lower box produces some opening) they start flying out which absolutely is not a case in the lab since soil in the lab doesn't come out from the lower box (after opening and after passing the thickness of upper shear box). I have the same issue when I try to compact soil layers in a container using UCM.

I think if I post a part of my script it might be more helpful. You won't be able to run it since you don't have the c++ codes for my contact model nor the packing I use for spheres but it gives some insights about what I am trying to do

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

from yade import ymport, utils, plot,export,qt
import time

################## parameters definition
Packing='/home/ngoudarz/Desktop/directShearSpheres7500Reload'
X=50e-3
Y=15e-3
Z=50e-3
S0=X*Y

dampingCoeff=0.5
dtCoeff=0.5
normalStress=75000
normalVel=0.1 # 0.001 for 100kPa // optimized for normalVEL=normalSTRESS/1e8?
shearVel=1*normalVel # try different values to ensure quasi-static conditions
intR=1.263

soilYoung = 1e8
steelYoung=210e9
soilPoisson=1.5
steelPoisson=.25
soilDensity = 2600
steelDensity = 0
depositionSoilFrictAngle = 0
experimentSoilFrictionAngle=0.5
steelFrictAngle = 0.0
soilBeta = 0.0
soilXic = 0.0
iterMax=400000

sphereMat = O.materials.append(FrictMat(young=soilYoung,poisson=soilPoisson,frictionAngle=0.523599, density=soilDensity,label='sphereMat'))
wallMat = O.materials.append(FrictMat(young=steelYoung,poisson=steelPoisson,frictionAngle=0, density=steelDensity,label='wallMat'))

for m in O.materials:
 print m.id
O.bodies.append(geom.facetBox((25e-3,7.5e-3,(5e-3)),(75e-3,7.5e-3,20e-3),wallMask=8,material=wallMat,wire=True,color=(0.2,0.5,1)))
lowerBackPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(5e-3)),(75e-3,7.5e-3,20e-3),wallMask=4,material=wallMat,wire=True,color=(0.2,0.5,1)))
lowerFrontPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(45.e-3)),(75e-3,7.5e-3,20e-3),wallMask=8,material=wallMat,wire=True,color=(0.2,0.5,1)))
upperBackPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(45.e-3)),(75e-3,7.5e-3,20e-3),wallMask=4,material=wallMat,wire=True,color=(0.2,0.5,1)))
upperFrontPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(5e-3-1.5e-3)),(25e-3,22.5e-3,20e-3),wallMask=2,material=wallMat,wire=True,color=(0.2,0.5,1)))
lowerRightPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(5e-3-1.5e-3)),(25e-3,22.5e-3,20e-3),wallMask=1,material=wallMat,wire=True,color=(0.2,0.5,1)))
lowerLeftPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(45.e-3+1.5e-3)),(25e-3,22.5e-3,20e-3),wallMask=2,material=wallMat,wire=True,color=(0.2,0.5,1)))
upperRightPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,(45.e-3+1.5e-3)),(25e-3,22.5e-3,20e-3),wallMask=1,material=wallMat,wire=True,color=(0.2,0.5,1)))
upperLeftPlane=O.bodies[-1]

O.bodies.append(geom.facetBox((-25e-3,7.5e-3,(12.5e-3+1.5e-3)),(25e-3,7.5e-3,12.5e-3),wallMask=32,material=wallMat,wire=False,color=(0.2,0.5,1)))
leftExtensionPlaneTop=O.bodies[-1]

O.bodies.append(geom.facetBox((-25e-3,7.5e-3,(12.5e-3-1.5e-3)),(25e-3,7.5e-3,12.5e-3),wallMask=32,material=wallMat,wire=False,color=(0.2,0.5,1)))
leftExtensionPlaneBottom=O.bodies[-1]

O.bodies.append(geom.facetBox((75e-3,7.5e-3,(12.5e-3-1.5e-3)),(25e-3,7.5e-3,12.5e-3),wallMask=32,material=wallMat,wire=False,color=(0.2,0.5,1)))
RightExtensionPlaneBottom=O.bodies[-1]

O.bodies.append(geom.facetBox((75e-3,7.5e-3,(12.5e-3+1.5e-3)),(25e-3,7.5e-3,12.5e-3),wallMask=32,material=wallMat,wire=False,color=(0.2,0.5,1)))
RightExtensionPlaneTop=O.bodies[-1]

O.bodies.append(geom.facetBox((25e-3,7.5e-3,25e-3),(75e-3,22.5e-3,25e-3),wallMask=16,material=wallMat,wire=True,color=(0.2,0.5,1)))
bottomPlate=O.bodies[-1]

id_spheres=O.bodies.append(ymport.text(Packing+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat,color=(0.75,0.1,0)))

for b in O.bodies:
   if isinstance (b.shape,Sphere):
      if (b.state.pos[2]-b.shape.radius)>55e-3:
         O.bodies.erase(b.id)

#cutMax_Z=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])

O.bodies.append(geom.facetBox((25e-3,7.5e-3,58e-3/2),(75e-3,22.5e-3,58e-3/2),material=wallMat,wallMask=32,wire=False,color=(1,1,1)))
topPlate=O.bodies[-1]

O.bodies.append(geom.facetBox((75e-3,7.5e-3,25e-3),(25e-3,7.5e-3,1.5e-3),material=wallMat,wallMask=1,wire=False,color=(1,1,1)))
rightGapPlate=O.bodies[-1]

O.bodies.append(geom.facetBox((-25e-3,7.5e-3,25e-3),(25e-3,7.5e-3,1.5e-3),material=wallMat,wallMask=2,wire=False,color=(1,1,1)))
leftGapPlate=O.bodies[-1]
#O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
#global topPlate
#topPlate=O.bodies[-1]
facet8LowerRight=O.bodies[8]
facet9LowerRight=O.bodies[9]
facet10LowerLeft=O.bodies[10]
facet11LowerLeft=O.bodies[11]
facet14UpperLeft=O.bodies[14]
facet15UpperLeft=O.bodies[15]
facet18LeftExtensionBottom=O.bodies[18]
facet19LeftExtensionBottom=O.bodies[19]
facet16LeftExtensionTop=O.bodies[16]
facet17LeftExtensionTop=O.bodies[17]
facet20RightExtensionBottom=O.bodies[20]
facet21RightExtensionBottom=O.bodies[21]
facet22RightExtensionTop=O.bodies[22]
facet23RightExtensionTop=O.bodies[23]
facet7526TopPlate=O.bodies[7526]
facet7527TopPlate=O.bodies[7527]

O.engines=[
   ForceResetter(),
   # sphere, facet
   InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='aabb'),Bo1_Facet_Aabb()],verletDist=(2.0e-3),label='collider',ompThreads=1),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='ss2d3dg'),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_JiangPhys(BETA=soilBeta,xIc=soilXic,label='sphereToSphereInteractionPhys'),Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=0,label='sphereToContainerInteractionPhys')],
      [Law2_ScGeom_JiangPhys_Jiang(includeRollResistMoment=False,includeTwistResistMoment=False,label='spherInteractionLaw'),Law2_ScGeom_FrictPhys_CundallStrack(label='sphereToContainerInteractionLaw')]
    ),
   NewtonIntegrator(gravity=(0.,0.,0.),damping=0.5,label='damper'),
   # the label creates an automatic variable referring to this engine
   # we use it below to change its attributes from the functions called
   GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),timestepSafetyCoefficient=dtCoeff),
   TranslationEngine(ids=[7526,7527],translationAxis=(0.,0.,-1.0),velocity=0,label='zTranslation'),
   TranslationEngine(ids=[8,9,10,11,18,19,20,21,7530,7531],translationAxis=(-1.0,0.,0.),velocity=0,label='xTranslation'),
   PyRunner(iterPeriod=1,initRun=True,command='normalCompaction()',label='checker'),
   NewtonIntegrator(damping=dampingCoeff,gravity=(0,0,0),label='damper'),
   PyRunner(iterPeriod=iterMax/1000,initRun=True,command='dataCollector()'),
]

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

def normalCompaction():
  global px0,pz0,sigmaN,n,Fn1,facet14UpperLeft,facet15UpperLeft,facet8LowerRight,facet9LowerRight,facet10LowerLeft,facet11LowerLeft,Fs11Norm,Fs12Norm,Fs1Norm
  Fn11=abs(O.forces.f(facet7526TopPlate.id)[2])
  Fn12=abs(O.forces.f(facet7527TopPlate.id)[2])
  Fn1=Fn11+Fn12
  sigmaN=Fn1/(S0) # necessary?
  Fs11Norm=abs(O.forces.f(facet14UpperLeft.id)[0])
  Fs12Norm=abs(O.forces.f(facet15UpperLeft.id)[0])
  Fs1Norm=Fs11Norm+Fs12Norm
  print "\r Fn1: ",Fn1,"Fn11: ",Fn11,"Fn12: ",Fn12,"sigmaN: ",sigmaN,"normalStress: ",normalStress,"zTranslation.velocity: ",zTranslation.velocity,"unbalanced: ",unbalancedForce()
  if zTranslation.velocity<normalVel:
     zTranslation.velocity+=normalVel/1000
  if sigmaN>(0.9*normalStress):
     zTranslation.velocity=normalVel*((normalStress-sigmaN)/normalStress)
  if abs((normalStress-sigmaN)/normalStress)<0.001 and unbalancedForce()<0.05:
     zTranslation.velocity=0.
     checker.command='directShear()'
       #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,0,1))/S0

### add engine and initialisation
px0=facet8LowerRight.state.pos[0]
pz0=facet7526TopPlate.state.pos[2]
def directShear():
  damper.gravity=(0,0,0)
  print 'shearing! || iteration=', O.iter
  #zTranslation.velocity=0. # if you want a constant normal diaplacement control
  zTranslation.velocity=50*normalVel*((normalStress-sigmaN)/normalStress) # not good enough because not general (coeff needs to be adjusted)
  if xTranslation.velocity<shearVel:
     xTranslation.velocity+=(shearVel/1000)
  if abs(facet8LowerRight.state.pos[0]-px0)>=(0.15*X):
     xTranslation.velocity=0.
     O.pause()

def dataCollector():
  global Xdispl,Zdispl,tauTot1,tauTot2,Fs1Net,Fs2,Fs11Tot,Fs12Tot,Fs1Tot
  Fs11Tot=abs(O.forces.f(facet14UpperLeft.id)[0])
  Fs12Tot=abs(O.forces.f(facet15UpperLeft.id)[0])
  Fs1Tot=Fs11Tot+Fs12Tot
  Fs1Net=abs(Fs1Tot-Fs1Norm)
  #Fs2=abs(O.forces.f(lowerRightPlane.id)[0])
  Xdispl=abs(facet8LowerRight.state.pos[0]-px0)
  Zdispl=abs(facet7526TopPlate.state.pos[2]-pz0)
  XRemained=abs(facet14UpperLeft.state.pos[0]-facet8LowerRight.state.pos[0])
  tauTot1=Fs1Net/(XRemained*Y)
  tauTot2=Fs1Net/(S0)
  print "\r Fs11Tot: ",Fs11Tot,"Fs12Tot: ",Fs12Tot,"Fs1Tot: ",Fs1Tot,"Fs1Net: ",Fs1Net,
  yade.plot.addData({'t':O.time,'i':O.iter,'i_':O.iter,'i__':O.iter,'Xdispl':Xdispl,'Xdispl_':Xdispl,'Xdispl__':Xdispl,'sigmaN':sigmaN,'Fn1':Fn1,'tauTot1':tauTot1,'tauTot2':tauTot2,'Fs1Norm':Fs1Norm,'Fs1Tot':Fs1Tot,'Fs1Net':Fs1Net,'Zdispl':Zdispl,'unbF':utils.unbalancedForce(),'Ek':utils.kineticEnergy()})
from yade import plot
plot.plots={'i':('sigmaN',),'Xdispl':('Fs1Net',),'i__':('Zdispl',),'Xdispl_':('tauTot1',),'Xdispl__':('tauTot2',),}
plot.plot()

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

Honestly I did not go in details through your whole post :-) but I'd still make the following comments

1- I think you could consider that in YADE / DEM, properties for a single entity (a Discrete Element) do not exist outside mass and inertia. Then, the question is more about : are the soil-container contact properties really different from the soil-soil ones ?

And if yes, does it really matter to define one external layer of interactions (the soil-container one) that has very different properties than all other interactions in the bulk ?

2- Checking the shape of interacting bodies can be done in the source code of one single Ip2 if you need to. You got one example at [*]

3- Are you then interested in shear stress (original post) or normal force and normal stress (post #2) ?

4- I still think material flowing outside a experimental or numerical device is a problem :-) If I understood correctly, I would just extend my plates if I were you (you may give a look to SimpleShear FileGenerator from YADE graphical user interface)

Jérôme

[*] https://github.com/yade/trunk/blob/4167bd08a304167375ac38bca441aa30e8d4d673/pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp#L101 to L103