JCFpmMat parameters define

Asked by Aaron Liu on 2020-11-29

Dear Yade users,

I found the JCFpmMat and related functors might be good to simulate bonded particle problems(https://answers.launchpad.net/yade/+question/635871). I was trying to simulate a uniaxial compression of a cubic brick that is making by using soft bonding material connecting hard sand particles(the mechanical parameters could be given). So my brick is made of two parts: soft bonding material and hard particles. Here are my questions.

In JCFpmMat (https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.JCFpmMat), it defines a lot of stiffness, strength parameters which are really confusing.
1. What are the differences between cohesion and jointCohesion?
2. What is the meaning of joint surface, does it describe the strength/stiffness of debonding? if so, why do we have to define the stiffness of debonding? why not just define the stiffness of the bonding material itself?
3. What is the young in my case, it depends on IP2 functor, but I did not find related information. Does it mean the Young modulus of particles?

from __future__ import division
from __future__ import print_function

from future import standard_library
standard_library.install_aliases()
from yade import plot,pack,timing
import time, sys, os, copy
readParamsFromTable(noTableOk=True, # unknownOk=True,
 young=1e6,
 poisson=.6,

 sigmaT=3.5e6,
 frictionAngle=radians(20),
 epsCrackOnset=1e-4,
 relDuctility=30,

 intRadius=1.5,
 dtSafety=.8,
 damping=0.4,
 strainRateTension=-1,
 strainRateCompression=-.5,
 setSpeeds=True,

 specimenLength=.04,
 sphereRadius=0.0005,

 # isotropic confinement (should be negative)
 isoPrestress=0,
 jointDil=radians(0),
 jointFrict=radians(20),
 OUT='compressionTest_JCFPM'

)

from yade.params.table import *
#alphaKr=30,alphaKtw=30,
#O.materials.append(JCFpmMat(type=1,young=1e10,frictionAngle=radians(30),density=2600,poisson=0.3,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1e7,jointShearStiffness=1e7,jointCohesion=1e6,jointFrictionAngle=atan(0.8),jointDilationAngle=radians(0))
O.materials.append(JCFpmMat(type=1,young=1e6,frictionAngle=radians(30),density=3000,poisson=0.3,tensileStrength=3e7,cohesion=1e6,jointNormalStiffness=1e7,jointShearStiffness=1e4,jointCohesion=1e4,jointFrictionAngle=jointFrict,jointDilationAngle=jointDil))

#O.materials.append(CohFrictMat(young=young,poisson=poisson,density=2600,frictionAngle=radians(30),normalCohesion=1e3,shearCohesion=1e3,momentRotationLaw=True,etaRoll=0.1,label='spheres'))
#concreteId=O.materials.append(FrictMat(young=15e6,poisson=.4,frictionAngle=0,density=0,label='frictionlessWalls'))

sp=pack.randomDensePack(pack.inAlignedBox((0,0,0),(specimenLength,specimenLength,specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation()
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=0.9*PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(label='is2aabb'),],verletDist=.05*sphereRadius),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom()],
 [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
 [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw'),Law2_ScGeom_FrictPhys_CundallStrack()]

 ),
 NewtonIntegrator(damping=damping,label='damper'),
 #CpmStateUpdater(realPeriod=.5),
UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
PyRunner(virtPeriod=1e-3/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 #PyRunner(realPeriod=4å,command='stopIfDamaged()',label='damageChecker'),
#
]

def addPlotData():
 yade.plot.addData({'t':O.time,'i':O.iter,'eps':-(strainer.strain),'sigma':-(strainer.avgStress+isoPrestress),})
print(strainer.strain)
plot.plots={'eps':('sigma',)}
plot.plot(subPlots=False)

O.saveTmp('initial');

O.timingEnabled=False

O.run()
 #InteractionLoop(
 #
# ),
#O.materials.append(JCFpmMat(type=1,young=1e8,frictionAngle=radians(30),density=3000,poisson=0.3,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1e7,jointShearStiffness=1e7,jointCohesion=1e6,jointFrictionAngle=atan(0.8),jointDilationAngle=radians(0))

Cheers,
Aaron

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2020-12-10
Last query:
2020-12-10
Last reply:
2020-11-30
Best Jérôme Duriez (jduriez) said : #1

Hi,

JCFpm parameters should be described e.g. in [Duriez2016] from https://www.yade-dem.org/doc/publications.html (with free copy at https://hal.archives-ouvertes.fr/hal-01865371), see § 3. for all "joint*" parameters, or 2.1 for the rest.

YADE "young" is denoted "Y" there, it is a particle-scale stiffness parameter, with wich the macroscopic Young modulus relates. It is an user input, what depends on Ip2 is its interpretation. (with your Ip2, interpretation is the one of the previous sentence)

Luc Scholtès (luc) said : #2

Hi Aaron,

In JCFpmMat (https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.JCFpmMat), it defines a lot of stiffness, strength parameters which are really confusing.

-> Well, YADE is a collaborative project so if you find things confusing, feel free to propose improvements. This being said (...), JCFPM was designed to deal with jointed materials, i.e., cohesive materials with pre-existing structures, like, for instance, fractured rocks (rocks with pre-existing structural discontinuities/fractures). Hence, all parameters are "doubled" so that the two constitutive component of a fractured rock (i.e., matrix and fractures) can have different mechanical properties. If the contact/bond between 2 interacting particles belong to a fracture plane, its properties are defined through the "jointDedicatedProperty". If the contact/bond between 2 interacting particles belong to the matrix, its properties are defined through the "dedicatedProperty". If your material does not contain any pre-existing structures, you can just forget about all the jointDedicatedProperties". Why did we bother with jointDedicatedProperties in JCFPM? Because we wanted to be able to reorientate the contact/bond plane accordingly to the discontinuity/fracture plane.

I strongly suggest you read [1].

1. What are the differences between cohesion and jointCohesion?

-> As mentioned above, "cohesion" is for interactions within the rock matrix. "jointCohesion" is for interactions along pre-existing discontinuities/fractures. The same logic applies to tensileStrength/jointTensileStrength, frictionAngle/jointFrictionAngle, YoungModulus/jointNormalStiffness, Poisson/jointShearStiffness. Then, you have extra parameters like jointDilation which enables to define the dilation of jointed contacts following the physics of rock joints.

2. What is the meaning of joint surface, does it describe the strength/stiffness of debonding? if so, why do we have to define the stiffness of debonding? why not just define the stiffness of the bonding material itself?

-> are you referring to the crossSection variable? I don't see any jointSurface variable? Or are you referring to the "jointNormal" variable? If that's the case, jointNormal is a vector that defines the orientation of the pre-exisitng discontinuity/fracture (if there is one)

3. What is the young in my case, it depends on IP2 functor, but I did not find related information. Does it mean the Young modulus of particles?

-> The name "Young" is confusing, I agree with you. It was defined at the very first stage of the code development and it kept going on afterward so we just have to deal with it... Actually, it directly relates to the macroscopic Young modulus. ANW, it defines the normal stiffness of contacts (https://yade-dem.org/doc/formulation.html?highlight=young#normal-stiffness). In JCFPM material, it defines the normal stiffness of contacts within the rock matrix (notOnJoint).

In your simulation, I'd suggest to use only one material type with 2 sets of properties insteaf of using both JCFpmMat and CohFrictMat.

Please find below a script that simulates a uniaxial compression tests on an intact matrial using JCFPM.

Cheers,

Luc

[1] https://www.sciencedirect.com/science/article/abs/pii/S1365160912000391: Luc Scholtès, Frédéric-Victor Donzé, Modelling progressive failure in fractured rock masses using a 3D discrete element method, International Journal of Rock Mechanics and Mining Sciences, Volume 52, 2012, Pages 18-30

----------

from yade import pack, plot

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

#### Name of output files
OUT='compressionTest'

#### Simulation Control
rate=-0.05 #deformation rate
iterMax=200000 # maximum number of iterations
saveVTK=int(iterMax/5.) # saving output files for paraview

#### Microproperties (interparticle parameters)
intR=1.4 # allows near neighbour interaction (can be adjusted for every packing / the bigger -> the more brittle / careful when intR is too large -> bonds can be created "over" particles) -> intR can be calibrated to reach a certain coordination number K (see calculation on line 115)
DENS=3000 # this one can be adjusted for different reasons (porosity of packing vs porosity of material / increase time step (no gravity -> no real effect on the result)
YOUNG=10e9 # this one controls the Young's modulus of the material
ALPHA=0.15 # this one controls the material Poisson's ratio of the material
TENS=3e6 # this one controls the tensile strength UTS of the material
COH=30e6 # this one controls the compressive strength UCS of the material, more precisely, the ratio UCS/UTS (from my experience: COH should be >= to TENS, >= 10*TENS for competent materials like concrete)
FRICT=30 # this one controls the slope of the failure envelope (effect mainly visible on triaxial compression tests)

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

#### create the specimen
L=0.10
D=0.05
pred=pack.inCylinder((0,0,0),(0,0,L),D/2.)
O.bodies.append(pack.regularHexa(pred,radius=D/20.,gap=0.,material=sphereMat)) # regular packings should be avoided as failure is controlled by particle arrangement

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

################# DEM loop + 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.8,defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
 VTKRecorder(iterPeriod=1,initRun=True,fileName=OUT+'-',recorders='spheres','jcfpm','cracks','bstresses'],Key=OUT,dead=1,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,
         'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

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

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();

### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

vtk.dead=0
O.step()

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

Aaron Liu (voyagening) said : #3

Hi Jérôme and Luc,

Thank u so much for the kind explanation. I will take a research on it further.

Best,
Aaron