yade-users team mailing list archive

[Luc Scholtès <question693282@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Needs information => Answered

Luc Scholtès proposed the following answer:
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))

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