Contact sliding and rolling friction values different than assigned?

Asked by Kevin Kuei

Hello,

I'm working off of a modified tutorial example. Code provided below for reference.

In summary:
- I am assigning a lower initial sliding and rolling friction during compaction, then bumping it up before shear.
    - Initial rolling/sliding friction = 0.0
    - Final rolling/sliding friction = radians(30) = 0.52359...
-Note the docs (help(O.interactions[index].phys)) say that maxRollPl is the coefficient of rollingfriction, so I am equating that essentially with etaRoll (maybe this is not exactly correct?).
- If I inspect the phys instances of O.interactions at the end of my simulation, I generally see that:
    - sliding friction = 0.577 = 33 deg
    - rolling friction < 0.577 overall, but sometimes as small as 4 deg
- I would expect them all to report 0.52359.., however, it is clearly different.

Just want to understand what I might be missing here.

Cheers,
Kevin

truncated output for contact friction at end of test:
```
for i in O.interactions:
   print(i.phys.maxRollPl, i.phys.tangensOfFrictionAngle ):

0.5773502691896257 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.04720909644502941 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.038848741799341516 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.04013287641764167 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.5773502691896257 0.5773502691896257
0.03819576331427891 0.5773502691896257
0.04293328164040254 0.5773502691896257
0.04918282798624137 0.5773502691896257
0.03809146798509421 0.5773502691896257
...
```

working example code:
```
from __future__ import print_function

# generate loose packing
from yade import pack, qt, plot

sigmaIso = -1e5
initialFric = 0
finalFric = 30

## define materials
O.materials.append(CohFrictMat(
        density=2650, # Density [kg/m3]
        young=1e7, # Particle modulus [Pa]
        poisson=.3, # Ks/Kn ratio
        frictionAngle=radians(initialFric), # Local friction [rad]
        isCohesive=False, # Turn off adhesion
        momentRotationLaw=True, # Turn on rotational stiffness
        etaRoll=radians(initialFric), # Rotational friction [rad]
        etaTwist=0, # Turn off twisting
        label="granr" # Material label
        ))

O.periodic = True
sp = pack.SpherePack()
## uniform distribution
sp.makeCloud((0, 0, 0), (4, 4, 4), rMean=.1, rRelFuzz=.3, periodic=True)

# setup periodic boundary, insert the packing
sp.toSimulation(material="granr")

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],
                        [Ip2_FrictMat_FrictMat_FrictPhys(),
                         Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                        [Law2_ScGeom_FrictPhys_CundallStrack(),
                         Law2_ScGeom6D_CohFrictPhys_CohesionMoment()],
                        label='interactionLoop'),
        GlobalStiffnessTimeStepper(
                # use adaptive stiffness-based timestepper
                timestepSafetyCoefficient=0.8,
                timeStepUpdateInterval=100,
                # set to True for domain decomp
                parallelMode=False,
                label='timeStepper'
        ),
        NewtonIntegrator(damping=.2),
        PeriTriaxController(
                label='triax',
                # specify target values and whether they are strains or stresses
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=7,
                # type of servo-control
                dynCell=True,
                maxStrainRate=(0.01, 0.01, 0.01),
                # wait until the unbalanced force goes below this value
                maxUnbalanced=1e-5,
                relStressTol=1e-5,
                # call this function when goal is reached and the packing is stable
                doneHook='compactionFinished()'
        ),
        PyRunner(command='addPlotData()', iterPeriod=100),
]

def addPlotData():
 plot.addData(
         unbalanced=unbalancedForce(),
         i=O.iter,
         sxx=triax.stress[0],
         syy=triax.stress[1],
         szz=triax.stress[2],
         exx=triax.strain[0],
         eyy=triax.strain[1],
         ezz=triax.strain[2],
         # save all available energy data
         Etot=O.energy.total(),
         **O.energy
 )

# enable energy tracking in the code
O.trackEnergy = True

# define what to plot
plot.plots = {
        'i': ('unbalanced',),
        'i ': ('sxx', 'syy', 'szz'),
        ' i': ('exx', 'eyy', 'ezz'),
        # energy plot
        ' i ': (O.energy.keys, None, 'Etot'),
}
# show the plot
plot.plot()

def compactionFinished():
    ## reassign friction values for shear
    # for future contacts change material
 O.materials[0].frictionAngle = radians(finalFric) # radians
 O.materials[0].etaRoll = radians(finalFric)
 # for existing contacts, set contact friction directly
 for i in O.interactions:
     i.phys.tangensOfFrictionAngle = tan(radians(finalFric))
     i.phys.maxRollPl = tan(radians(finalFric))

    # 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
 triax.goal = (sigmaIso, sigmaIso, -.2)
 triax.stressMask = 3
 # allow faster deformation along x,y to better maintain stresses
 triax.maxStrainRate = (1., 1., .01)
 # next time, call triaxFinished instead of compactionFinished
 triax.doneHook = 'triaxFinished()'
 # do not wait for stabilization before calling triaxFinished
 triax.maxUnbalanced = 10

def triaxFinished():
 print('Finished')
 O.pause()

```

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,

I probably did not get completely (all) your question(s) but I can tell that CohFrictPhys.maxRollPl is not supposed to be equal to CohFrictMat.etaRoll. The former has a length unit while the later is dimensionless.

While it may not be really explained in the doc, code-source might help:

https://gitlab.com/yade-dev/trunk/-/blob/91151dd317a97f5ac11a8fbeafb565c853899a67/pkg/dem/CohesiveFrictionalContactLaw.cpp#L321 (current YADE version)

Can you help with this problem?

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

To post a message you must log in.