Friction angle between synthetic rock and Cutter

Asked by kalogeropoulos on 2019-01-27

Hi to all,

I have written a code for the calculation of the cutting force, and the two others during a cutting test,

My question is : I am not sure i have put the friction angle between rock-cutter right.

As you will see i defined the two materials, and i put there the friction angle, one for sphere-sphere, and the other for facet-sphere. Is it right or i have to define it inside Ip2 or anywhere else??

Help!!!!!

### JCFPM script

# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from __future__ import division
from yade import pack, plot, geom, pack,utils
from yade import _packPredicates

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

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

#### Simulation Control
#rate=0.01#-0.01 #deformation rate
#iterMax=800000 # maximum number of iterations Doesnt Need anymore
saveVTK=4000 # saving output files for paraview

#### Material microproperties
intR=1.25 # allows near neighbour interaction (can be adjusted for every packing)
DENS=5290.26 #2674.3 toy petrwmatos## 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=42e9
FRICT=43.94
ALPHA=2.5
TENS=11e6
COH=1100e6

#### Cutter Properties
x_cutter = 6.35e-3 # Width along x-axis
y_cutter = 1e-3 # Height along y-axis
z_cutter = 5e-3 # Length along z-axis
DOC = 5e-3 # Depth of Cut

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
sample = O.materials.append(sphereMat())
#### create the specimen
X = 54.85e-3
Y=54.85e-3
Z=146.5e-3
pred= pack.inAlignedBox((0,0,0),(X,Y,Z))
sps=SpherePack()
sp=pack.randomDensePack(pred,spheresInCell=10000,radius=1.3725e-3,rRelFuzz=.34,memoizeDb='/tmp/07012019.sqlite',returnSpherePack=True)
sp.toSimulation(material=sample) # up to you to define another sample here, e.g., with randomDensePack or anything else.

#### create the cutter
O.materials.append(JCFpmMat(young=600e9,poisson=0.225,frictionAngle=radians(30),density=1470000,label='bx'))
bx = geom.facetBox(center=(27.425e-3,Y+5e-3,-5e-3),extents=(x_cutter,-(5e-3+DOC),z_cutter),orientation=Quaternion((1,0,0),pi/36),wallMask=8+32)
O.bodies.append(bx)

for facet in bx:
    facet.state.blockedDOFs='xyzXYZ'
    facet.state.vel=(0.0,0.,5.0) # velocity along x-axis
      #block e.w na metakineitai mono kata

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

################# ENGINES DEFINED HERE
## damping = 0.4

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb'),Bo1_Facet_Aabb()],verletDist=.05*1.3725e-3),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.6,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
     # VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks','facets'],Key=OUT,label='vtk'),
 PyRunner(iterPeriod=int(100), command='Stop()',label='data')
]

################# RECORDER DEFINED HERE
def recorder():
   global Fz, Fx, Fy
   global dz
   Fz = 0.0
   Fz = abs(sum(O.forces.f(facet.id)[2] for facet in bx))
   Fx = abs(sum(O.forces.f(facet.id)[0] for facet in bx))
   Fy = abs(sum(O.forces.f(facet.id)[1] for facet in bx))
   yade.plot.addData({'i':O.iter,
        'Fx':Fx,
        'Fy':Fy,
        'Fz':Fz,
        'dz': abs(bx[2].state.pos[2]),
        'tc':interactionLaw.nbTensCracks,
        'sc':interactionLaw.nbShearCracks,
        'te':interactionLaw.totalTensCracksE,
        'se':interactionLaw.totalShearCracksE,
        'unbF':utils.unbalancedForce()})
   plot.saveDataTxt(OUT)

def Stop():
    if bx[2].state.pos[2] >= Z+1e-3:
        O.pause()

# if you want to plot during simulation
plot.plots={'dz':('Fz')}
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
print 40*'-'
################# SIMULATION REALLY STARTS HERE
#O.run()

Thanks a Lot!!

A.K

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2019-01-28
Last query:
2019-01-28
Last reply:
2019-01-28
Jérôme Duriez (jduriez) said : #1

Hi,

Your script is way too long for me to look into it.
You can check yourself the friction angle you got at the interaction level, though.

Eg O.interactions[idSphere1,idSphere2].phys.tanFrictionAngle to compare with
O.interactions[idSphere1,idFacet].phys.tanFrictionAngle

You can also use "Inspect" from the Graphical User Interface.

Jérôme

kalogeropoulos (antoniskal) said : #2

In the above code, i don't know where to set the friction angle of the interaction cutter-rock.

Because Ip2_JCFpmMat_JCFpmMat_JCFpmPhys() doesn't support tanfrictionAngle, do i have to put and the

--> Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=30) in the engines???

Best Jérôme Duriez (jduriez) said : #3

No, JCF* model follows the same logic as Frict* : mechanical properties are stored/defined by the user in bodies' Material (JCFpmMat, or FrictMat, or ...), and just computed by the corresponding Ip2. (It's true the doc could explain this better ?)

What you should do here is to define two JCFpmMat : one for rock particles, one for cutter particles, differing in terms of friction angles (if that's your goal).

The friction angle of an interaction between two bodies with these two different JCFpmMat will be the minimum of the two friction angles [*], as computed by the Ip2 [*]

I would advice to start with just two particles.

[*] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/JointedCohesiveFrictionalPM.cpp#L537

kalogeropoulos (antoniskal) said : #4

Thanks Jérôme Duriez, that solved my question.