Perform UCS-TS test

Asked by kalogeropoulos on 2018-12-06

Hi to all,

currently i am trying to perform simulation of UCS and TS as in the paper <<A DEM model for soft and hard rocks: Role of grain interlocking on strength >> of Luc Scholtes and Frédéric-Victor Donzé.

I took the file from trunk/examples/concrete/uniax.py and change the model from CpmMat to JCFpmMat.
Also i took the microparameters from table 2 in the above paper and put it in my script.

Some Remarks i saw in my simulation:
As i put intRadius=1 as in table 2, my Coordination number is 3.37 byt in paper is 6.1
Number of elements: 17412
And Fc = 1.36e06 .
Am i doing something wrong?
 Please have a look in my code :

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing
import time, sys, os, copy

#import matplotlib
#matplotlib.rc('text',usetex=True)

# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
    density=2640,
 young=230e9,
 poisson=2.5,
 frictionAngle=25,
 tensileStrength=22e6,
 cohesion=220e6, #10*tensileStrength,

 sphereRadius=1.25e-3,
 rRelFuzz=0.34,

 intRadius=1,
 dtSafety=.8,
 damping=0.6,
 strainRateTension=.05,
 strainRateCompression=.5,
 setSpeeds=True,
 # 1=tension, 2=compression (ANDed; 3=both)
 doModes=2,

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']

sample = O.materials.append(JCFpmMat(young=young, poisson=poisson, frictionAngle=radians(frictionAngle), cohesion=cohesion, tensileStrength=tensileStrength,density = density, label='spheres'))

sps=SpherePack()

sp=pack.randomDensePack(pack.inAlignedBox((-25e-3,-25e-3,-62.5e-3),(25e-3,25e-3,62.5e-3)),spheresInCell=2000,radius=sphereRadius,rRelFuzz=rRelFuzz,memoizeDb=False,returnSpherePack=True) #'/tmp/triaxPackCache.sqlite'

sp.toSimulation(color=(1.09,1.0,0.0),material=sample) #color=(1.09,1.0,0.0)

bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt

mm,mx=[pt[axis] for pt in aabbExtrema()]
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb')],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
        [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(Key='cracks',recordCracks=True,cracksFileExist=True,label='interactionLaw')],
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
    NewtonIntegrator(damping=damping,label='damper'),
        VTKRecorder(fileName='3d-vtk-',recorders=['intr','cracks','jcfpm','facets','spheres','colors'],Key='cracks', realPeriod=50),
 PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]

plot.plots={'eps':('sigma',)}

O.saveTmp('initial');

O.timingEnabled=False

global mode
mode='tension' if doModes & 1 else 'compression'

def initTest():
 global mode
 print "init"
 if O.iter>0:
  O.wait();
  O.loadTmp('initial')
  print "Reversing plot data"; plot.reverseData()
 else: plot.plot(subPlots=False)
 strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
 try:
  from yade import qt
  renderer=qt.Renderer()
  renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
 except ImportError: pass
 print "init done, will now run."
 O.step(); # to create initial contacts
 # now reset the interaction radius and go ahead
 ss2sc.interactionDetectionFactor=1.
 is2aabb.aabbEnlargeFactor=1.

 O.run()

def stopIfDamaged():
 global mode
 if O.iter<2 or not plot.data.has_key('sigma'): return # do nothing at the very beginning
 sigma,eps=plot.data['sigma'],plot.data['eps']
 extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
 minMaxRatio=0.5 if mode=='tension' else 0.5
 if extremum==0: return
 import sys; sys.stdout.flush()
 if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
  if mode=='tension' and doModes & 2: # only if compression is enabled
   mode='compression'
   O.save('uniax-tension.yade.gz')
   print "Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)"
   print "Damaged, switching to compression... "; O.pause()
   # important! initTest must be launched in a separate thread;
   # otherwise O.load would wait for the iteration to finish,
   # but it would wait for initTest to return and deadlock would result
   import thread; thread.start_new_thread(initTest,())
   return
  else:
   print "Damaged, stopping."
   ft,fc=max(sigma),min(sigma)
   print 'Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft))
   title=O.tags['description'] if 'description' in O.tags.keys() else O.tags['params']
   print 'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
   print 'Bye.'
   O.pause()
   #sys.exit(0) # results in some threading exception

def addPlotData():
 yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
  'sigma.25':forcesOnCoordPlane(coord_25,axis)[axis]/area_25+isoPrestress,
  'sigma.50':forcesOnCoordPlane(coord_50,axis)[axis]/area_50+isoPrestress,
  'sigma.75':forcesOnCoordPlane(coord_75,axis)[axis]/area_75+isoPrestress,
  })
plot.plot(subPlots=False)
#O.run()
initTest()

Thank You!!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Luc Scholtès
Solved:
2018-12-07
Last query:
2018-12-07
Last reply:
2018-12-06
Best Luc Scholtès (luc) said : #1

Hi,

For your intention, please find below a script that should do what you want.

Regarding your point about the coordination number, you have to keep in mind that it is directly related to the sample. For instance, a loose sample will present a smaller coordination number than a dense sample if you define the same interaction range in both cases. In the paper you mentioned, I generated my samples beforehand by compressing (hydrostatically) clouds of polydisperse particles. I choose the set of parameters (contact stiffnesses and confining pressure) so that the compacity of the sample was "high enough" and the interpenetration between particles "limited". Up to you to choose a generation procedure that suits your needs. IMO, for rock like materials, dense samples give more consistent results. You may need to experience different generation procedures to make your choice. randomDensePack is an option but I never found it satisfactory.

Hope it helps

Luc

### JCFPM script

# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import pack, plot

################# SIMULATIONS DEFINED HERE

#### packing (previously constructed)
OUT='compressionTest_JCFPM'

#### Simulation Control
rate=-0.01 #deformation rate
iterMax=10000 # maximum number of iterations
saveVTK=2000 # saving output files for paraview

#### Material microproperties
intR=1.1 # allows near neighbour interaction (can be adjusted for every packing)
DENS=2500 # could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing porosity as to be computed?
YOUNG=20e9
FRICT=7
ALPHA=0.1
TENS=1e6
COH=1e6

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)

#### create the specimen
pred=pack.inCylinder((0,0,0),(0,1,0),0.25)
O.bodies.append(pack.regularHexa(pred,radius=0.025,gap=0.,material=sphereMat)) # up to you to define another sample here, e.g., with randomDensePack or anything else.

R=0
Rmax=0
nbSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/nbSpheres

print 'nbSpheres=',nbSpheres,' | Rmean=',Rmean

#### boundary condition (see utils.uniaxialTestFeatures
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
         'eps':strainer.strain,
         'sigma':strainer.avgStress,
         'tc':interactionLaw.nbTensCracks,
         'sc':interactionLaw.nbShearCracks,
         'te':interactionLaw.totalTensCracksE,
         'se':interactionLaw.totalShearCracksE,
         'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

# if you want to plot during simulation
plot.plots={'i':('sigma')}
plot.plot()

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

#### coordination number verification
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
    if i.phys.isCohesive :
      numCohesivelinks+=1

print "nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, " || Kcohesive=", 2.0*numCohesivelinks/nbSpheres

################# SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(iterMax)

kalogeropoulos (antoniskal) said : #2

Thank you Luc, that may solve my problem!!!

kalogeropoulos (antoniskal) said : #3

Thanks Luc Scholtès, that solved my question.