Does 'JCFpmPhys' have 'matchmaker'?

Asked by William

Hi, everyone!
I want to simulate the combination of granular particles and rocks. That is to say, one part of a sample is soil particles, and the other part is rock consist of the bonded particles.

Firstly, I import the particle using ‘ymport.text’. Then, I create two types of sphere materials (the tensileStrength and cohesion between the materials are different) and assign them to the particles in percentage terms.

Here is the mats:
mat1=JCFpmMat(type=1,young=Young,poisson=Poisson,frictionAngle=radians(compFricDegree),density=Density,tensileStrength=10e4,cohesion=45e4,label='spheres1')
mat_1=O.materials.append(mat1)
mat2=JCFpmMat(type=1,young=Young,poisson=Poisson,frictionAngle=radians(compFricDegree),density=Density,tensileStrength=0,cohesion=0,label='spheres2')
mat_2=O.materials.append(mat2)

How to set the tensileStrength and cohesion between the particles of different mats? For example, is the tensileStrength 10e4, 0 or other value between the particles of mat1 and mat2? Is the ‘matchmaker’ needed in the engine here?

Here is my script.

#### packing (previously constructed)
PACKING='Position.txt'

#### name of output files
OUT=PACKING+'_1MPa_r0.02'

from yade import pack, timing,ymport

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

## Sphere material
Young=300e6
Density=2650
Poisson=0.3
compFricDegree = 0
finalFricDegree = 27
AvgRadius=0.0002

intR=1
TENS=0.1e6
COH=0.45e6
TENS0=0
COH0=0

## wall material
WYoung=600e7
WPoisson=0.5
WDensity=0
WFrictionAngle=0

damp=0.7 # damping coefficient
stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below)

# corners of the initial packing
size=0.05
mn,mx=Vector3(0,0,0),Vector3(size,size,size)
#mn2,mx2=Vector3(0,0,0.0025),Vector3(size,1.5*size,0.0025)

O.materials.append(JCFpmMat(type=0,young = WYoung , poisson = WPoisson , frictionAngle = radians(WFrictionAngle) , density = WDensity,label='walls'))
#O.materials.append(JCFpmMat(type=1,young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree),density = Density,tensileStrength=TENS,cohesion=COH,label='spheres'))
walls=aabbWalls([mn,mx],thickness=0.001,material='walls')
wallIds=O.bodies.append(walls)

mat1=JCFpmMat(type=1,young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree),density = Density,tensileStrength=TENS,cohesion=COH,label='spheres1')
mat_1=O.materials.append(mat1)
mat2=JCFpmMat(type=1,young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree),density = Density,tensileStrength=TENS0,cohesion=COH0,label='spheres2')
mat_2=O.materials.append(mat2)

O.bodies.append(ymport.text(PACKING,scale=1.,shift=Vector3(0,0,0)))

sphere1=[]
sphere2=[]
for b in O.bodies:
  if not isinstance(b.shape,Sphere): # change material only on spheres
    continue
  if random.random() < 0.1:
    b.mat = mat1
    b.shape.color = (1,0,0)
    #b.state.mass*=mat1.density/mat1.density
    sphere1.append(b.id)
  else:
    b.mat = mat2
    b.shape.color = (0,1,1)
    #b.state.mass*=mat2.density/mat1.density
    sphere2.append(b.id)

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

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
    #wall_back_activated = False, #for 2d simulation
    #wall_front_activated = False,
    thickness = 0,
    stressMask = 7,
    internalCompaction=False, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

triax.goal1=triax.goal2=triax.goal3=-100000
setContactFriction(radians(finalFricDegree))
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
                [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
                [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
 #PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
]

while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print ('unbalanced force:',unb,' mean stress: ',triax.meanStress," porosity:",triax.porosity,)
  if unb<0.000001 and abs(-100000-triax.meanStress)/100000<0.001:
    break

O.save('2.yade.gz')
O.save('2.yade.bz2')

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jérôme Duriez (jduriez) said :
#1

Hi,

*Phys classes never (I think..) have MatchMaker instances in their attributes, they have actual interaction properties values, after computations by Ip2_*Phys. Maybe the corresponding doc, about Ip2_JCFpmMat_JCFpmMat_JCFpmPhys, is indeed imperfect regarding the mixture rule necessary for your case.

Source code (e.g., [*]) tells us the minimum Material value takes it all, and there will be no tensile strength between mat1- and mat2-material particles.

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

Can you help with this problem?

Provide an answer of your own, or ask William for more information if necessary.

To post a message you must log in.