# yade-users team mailing list archive

## [Question #695332]: Triaxial test with rolling resistance effect

Hello,

I am working on the behavior of granular material under the triaxial test. I want to find the effect of rolling resistance on the shear behavior of the material. I simulated as the instruction, but I have no idea why the rolling resistance does not affect the shear response. A very short version of my script is here. I was wondering if you could give me a suggestion if I am doing wrong in my code.

Best Regards,

Alireza

from datetime import datetime
import math
import pylab
from numpy import *

#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3, frictionAngle= 0.31, fragile=False, alphaKr=1.5, momentRotationLaw=True, etaRoll=0.5, label='aggregate-48'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1532.2, poisson=0.3, frictionAngle= 0.31, fragile=False, alphaKr=1.5, momentRotationLaw=True, etaRoll=0.5, label='aggregate-814'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=2e11, density=0, poisson=0.3, frictionAngle= 0, fragile=False, alphaKr=0, momentRotationLaw=True, etaRoll=0, label='wall-frictionless'))

#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(9e-2,9e-2,9e-2)],thickness=0.001,oversizeFactor=1.5,material='wall-frictionless')
wallIds=O.bodies.append(walls)

#===============================================================
#=================== define packing ============================
#===============================================================

nums=['t','t','t']

mats=['aggregate-48','aggregate-814','aggregate-814']

coke=((1.875e-3,500),(0.9475e-3,1074),(0.50125e-3,2090))

color=((0,0,1),(0,1,0),(1,0,0))

tolerance=[(0.4752e-3),(0.2055e-3),(0.1539e-3)]

for i in range(len(nums)):

nums[i]=pack.SpherePack()
nums[i].makeCloud((5e-3,5e-3,5e-3),(85e-3,85e-3,85e-3),rMean=coke[i][0],rRelFuzz=tolerance[i],num=coke[i][1])
O.bodies.append([utils.sphere(c,r,material=mats[i],color=color[i]) for c,r in nums[i]])

#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
),
TriaxialStateRecorder(iterPeriod=100000,file='WallStresses'),
VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=100000),
#        qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
NewtonIntegrator(damping=0.4,gravity=[0,0,-9.81]),
PyRunner(command='calm()',iterPeriod=10,label='calmEngine')
]

O.dt=6e-7

plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2], ev=(triax.strain[0]+triax.strain[1]+triax.strain[2]), sxx=-triax.stress(triax.wall_right_id)[0], syy=-triax.stress(triax.wall_top_id)[1], szz=-triax.stress(triax.wall_front_id)[2], q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]), p=triax.meanStress, R=triax.stress(triax.wall_front_id)[2]/sigmaIso, por=triax.porosity, i=O.iter,**O.energy,total=O.energy.total()),

#=================================================================
#=================  APPLYING CONFINING PRESSURE   ================
#=================================================================
sigmaIso=-1e5
triax=TriaxialStressController(
maxMultiplier=1.000,
finalMaxMultiplier=1.000,
thickness = 0,
internalCompaction=False,
max_vel=0.05,
)

O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True

while 1:
O.run(100, True)
unb=unbalancedForce()
meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
print('mean Stress:',triax.meanStress,'porosity:', triax.porosity,'Unbalancing:',unb)
if abs(triax.meanStress)<100:
print("###      calm engine is dead     ###")
if  triax.porosity<0.5 and abs(sigmaIso-triax.meanStress)/abs(sigmaIso)<0.01:
print('Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:', -triax.strain[2])
break
print("###      Isotropic state saved      ###")

#====================================================================
#====================   Piont A localization   ======================
#====================================================================

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width

#====================================================================
#====================================================================

triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = -0.05
triax.max_vel=0.05

O.trackEnergy=True
O.saveTmp()

#====================================================================
#==================     Micro strain   =====================
#====================================================================
TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.volume(10)
TW.setState(0)
step=0
while 1:
step +=1
O.run(100000,True)
print('Dev strain1:',-triax.strain[0], 'Dev strain 2:',-triax.strain[1], 'Dev strain 3:', -triax.strain[2])
TW.setState(1)
TW.defToVtk("strain_"+str(step)+".vtk")
if abs(triax.strain[2]) > 0.5:
#		makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
print(O.iter)
print('time is ', O.time)
print('porosity:',triax.porosity)
print('strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2])