yade-users team mailing list archive

[Kevin Kuei <question707239@xxxxxxxxxxxxxxxxxxxxx>]

New question #707239 on Yade:
https://answers.launchpad.net/yade/+question/707239

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()

```





-- 
You received this question notification because your team yade-users is
an answer contact for Yade.