Polyhedra splitting with Mohr-Coulomb-Weibull criterion

Dear Yade Users,

I am trying to understand and use an implementation of a combined Mohr-Coulomb–Weibull (MCW) failure criterion [1].

I have created an example based on the "splitter.py" example prepared by François Kneib. My example consists of a cube compressed between two walls. I created a polyhedraMat, which is splittable and has defined all the properties of Mohr-Coulomb–Weibull failure criterion. I set all the strengths to 1, probability P=0.63, and high Weibull modulus m=1000 (just in case). I expect to get strength close to 1 Pa. Instead of that, the strength is around 1.6 e5. It is independent of any parameter dedicated for MCW criterion. On the other hand, it depends on "strength" material property (it is not included in the example, but the default value is 100).

Yade version is 2018.02b

The major issue is:
- Do I have to activate the MCW criterion somehow (or deactivate the previous one [2])?

And there are two minor issues that I do not understand:
- Do I have to set "Wei_S0". From what I understood from [1] it is some general value that could be used if the tension and compression strengths were the same. But if I declare both strengths (sigmaCD, sigmaCZ), I expect that Wei_S0 is calculated based on them. Am I right?
- Since the default "strength" is 100, why do I get 1.6e5 in my "experiment"? The scale effect should even decrease strength. In the experiment, the contact force between the wall and the cube is divided by crosssection area (1m2) to calculate stress.

Cheers,
Karol

[1] Gladkyy, A., & Kuna, M. (2017). DEM simulation of polyhedral particle cracking using a combined Mohr-Coulomb–Weibull failure criterion. Granular Matter, 19(3), 41.
[2] Eliáš, J. (2014). Simulation of railway ballast using crushable polyhedral particles. Powder Technology, 264, 458-465.
---------------------------------------
from yade import polyhedra_utils
from yade import plot

# FUNCTION DENINITIONS
def myAddData():
 force = O.forces.f(1)[1]# cube wall interaction
 sig = force/(d**2)
 plot.addData(t = O.time, force = force, sig = sig)

d = 1 # cube length

#cube material
gravel1 = PolyhedraMat()
gravel1.IsSplitable = True
gravel1.Wei_P = 0.63
gravel1.Wei_m = 1000
gravel1.Wei_S0 = 1
gravel1.Wei_V0 = d**3
gravel1.sigmaCD = 1
gravel1.sigmaCZ = 1

#wall material
steel = PolyhedraMat()
steel.young = 1e10
steel.frictionAngle = 0

# create cube
p1 = polyhedra_utils.polyhedra(gravel1, v= [(0,0,0),(0,d,0),(d,d,0),(d,0,0),(0,0,d),(0,d,d),(d,d,d),(d,0,d)])

w1 = utils.wall(+d*1.001, axis=1, sense=-1, material=steel)
w2 = utils.wall(-0.001*d, axis=1, sense=+1, material=steel)

v = 5e-4
w1.state.vel = (0,-v,0)
O.bodies.append((p1,w1,w2))

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
   ),
   NewtonIntegrator(),
 PolyhedraSplitter(iterPeriod=100),
 PyRunner(command='myAddData()',iterPeriod=2000),
]

O.dt= 0.1*polyhedra_utils.PWaveTimeStep()

plot.plots={'t':'sig'}
plot.plot(subPlots =False)

try:
 from yade import qt
 qt.Controller()
 v = qt.View()
 v.ortho = True
except:
 pass

O.run(50000),