yade-users team mailing list archive

[Nima Goudarzi <question668213@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Nima Goudarzi is still having a problem:
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)

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.