[Question #707253]: Biaxial test with rolling friction

```New question #707253 on Yade:

Hello everyone,
I want to conduct a biaxial test for sands with rolling resistance (by CohFrictMat). For a dense packing, I need to set the rolling resistance and friction angle to 0 in the compaction period, and then activate them to the real value in shearing period.
The modification of friction angle can be achieved by "setContactFriction". But for rolling resistance, I'm not sure the following is the right approach: set the value of "alaways_use_moment_law" as "false" in compaction period and modify it as "true" in shearing period. This approach draws on the example named "triax-cohesive" [1]. And this method is not quite the same as the previous approaches in the Q&A (for example, [2]). Is there any difference between these two approaches?

####Here is my code####
from __future__ import print_function
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

O.materials.append(CohFrictMat(young=1.e9,poisson=.8,alphaKr=1.0,frictionAngle=.0,momentRotationLaw=True,etaRoll=0.3,label='balls'))

sp = pack.SpherePack()

dmax=0.030
size=5*dmax
thickness=size/3.0

psdSizes=[0.0020,0.0048,0.0076,0.0104,0.0132,0.0160,0.0188,0.0216,0.0244,0.0272,0.0300]
psdCumm =[0.0000,0.3228,0.4924,0.6081,0.6961,0.7671,0.8266,0.8778,0.9228,0.9629,1.0000]

sp.makeCloud(minCorner=(0,0,thickness/2),maxCorner=(size,size,thickness/2),psdSizes=psdSizes,psdCumm=psdCumm,num=1400,periodic=True,seed=1,distributeMass=True)
sp.toSimulation()
O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,thickness)    # RVE scale up
print(len(O.bodies))
for p in O.bodies:
p.state.blockedDOFs = 'zXY'
p.state.mass = 2650 * thickness * pi * p.shape.radius**2 # 0.1 = thickness of cylindrical particle   # RVE scale up
inertia = 0.5 * p.state.mass * p.shape.radius**2
p.state.inertia = (.5*inertia,.5*inertia,inertia)

O.dt = utils.PWaveTimeStep()
print(O.dt)

O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=False, useIncrementalForm=True)]
),
PeriTriaxController(
dynCell=True,
goal=(-4.e4,-4.e4,0),
relStressTol=.001,
maxUnbalanced=.001,
maxStrainRate=(.5,.5,.0),
doneHook='compactionFinished()',
label='biax'
),
NewtonIntegrator(damping=.1),
]

sx=biax.stress[0],
sy=biax.stress[1],
sz=biax.stress[2],
ex=biax.strain[0],
ey=biax.strain[1],
ez=biax.strain[2],
ev=biax.strain[0]+biax.strain[1]+biax.strain[2],
eyp=-biax.strain[1],
sxp=-biax.stress[0],
syp=-biax.stress[1],
i=O.iter,
unbalanced=utils.unbalancedForce(),
totalEnergy=O.energy.total(),
**O.energy  # plot all energies
)
# enable energy tracking in the code
O.trackEnergy=True

plot.plots = {
'i': (('unbalanced', 'go'), None, 'kinetic'),
' i': ('ex', 'ey', 'ez', None, 'sx', 'sy', 'sz'),
'eyp ': ('ev'),
' eyp ': ('sxp', 'syp', None, 'ev'),###'a', ' a', 'a ', ' a ' means the different position of each plot
}
plot.plot()

def compactionFinished():
print(biax.stress)
print(O.cell.hSize)
print(vr)
O.engines[2].lawDispatcher.functors[0].always_use_moment_law=True
# set the current cell configuration to be the reference one
O.cell.trsf=Matrix3.Identity
# change control type: keep constant confinement in x,y, 20% compression in z
biax.goal=(-4.e4,-0.1,0)
# allow faster deformation along x,y to better maintain stresses      # RVE scale up, y value should be scale down a lot
biax.maxStrainRate=(1.0,0.05,0.0)
# next time, call triaxFinished instead of compactionFinished
biax.doneHook='biaxFinished()'
# do not wait for stabilization before calling triaxFinished
biax.maxUnbalanced=10

def biaxFinished():
print(O.cell.hSize)
vr = yade._utils.voidratio2D(zlen=thickness)             # RVE scale up
print(vr)
print('Finished')
O.pause()