yade-users team mailing list archive

[Kevin Kuei <question708020@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello, 

I am having some difficulties with reassigning the rolling friction using the PeriTriaxController after compaction. 

In the example code below, I use an initial rolling friction 'eta_roll_init' of 0 (sliding fric also set to 0) to achieve a dense packing, then reassigning a new rolling friction 'eta_roll_final' at end of compaction. The issue seems to be that using different 'eta_roll_final' values yields identical specimen responses at the end of shearing. 

I am not sure exactly what the issue is and am  wondering if someone with experience with this model can comment if I've doing the material and existing contact re-assignment correctly?

For reference, I used these previous answers as a guide:
https://answers.launchpad.net/yade/+question/701184
https://answers.launchpad.net/yade/+question/251218
https://answers.launchpad.net/yade/+question/251218

My understanding is that for existing contacts at end of compaction, I need to only update/re-calculate the maxRollPl value which is calculated as:

min(rolling_friction * ball1_radius, rolling_friction * ball2_radius)


Example Code:

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

from __future__ import print_function
from yade import pack, qt, plot, export

# Define model parameters
density = 2650           # Density (kg/m3)
young = 1e8               # Contact modulus (Pa)
poisson = 0.25           # Shear to normal stiffness ratio, Ks/Kn
fric_final = 20.0         # Sliding friction (deg)
alpha_roll = 1.0         # Rolling resistance ratio, Kr/(R1*R2*Ks)
eta_roll_final = 1.0   # Rolling friction (rad)

# Define PSD and RVE parameters
d_min, d_max = 20*0.13/1e3, 20*0.3/1e3         # ..log-linear
Lx, Ly, Lz = 15*d_max, 15*d_max, 15*d_max
numbins = 25

# Define specimen preparation parameters
sigma_iso = -1e5
fric_init = 0.0
eta_roll_init = 0.0
packing = 0
max_strain_rate_iso = 0.20
max_unbalanced_iso = 1e-5
rel_stress_tol_iso = 1e-5

# Define controller parameters
strain_target = 0.20
max_strain_rate_horiz = 1.0
max_strain_rate_vertical = 0.20
max_unbalanced_force = 10


# Define materials
O.materials.append(CohFrictMat(
  density=density,                                  # Density [kg/m3]
  young=young,                                      # Particle modulus [Pa]
  poisson=poisson,                                # Ks/Kn ratio
  frictionAngle=radians(fric_init),    # Local friction [rad]
  isCohesive=False,                               # Turn off adhesion
  momentRotationLaw=True,           # Turn on rotational stiffness
  alphaKr=alpha_roll,                           # Dimensionless rolling stiffness
  alphaKtw=0,                                        # Dimensionless twist stiffness
  etaRoll=eta_roll_init,                       # Rotational friction [rad]
  etaTwist=0,                                         # Turn off twisting
  label="granr"                                     # Material label
  ))


# Define periodic boundaries
O.periodic = True

# Prepare specimen packing
sp = pack.SpherePack()
# ..log-linear distribution
start, end = numpy.log10(d_min), numpy.log10(d_max)
psdSizes = numpy.logspace(start, end, num=numbins, endpoint=True) 
psdCumm = numpy.linspace(0, 1., len(psdSizes))
sp.makeCloud((0, 0, 0), (Lx, Ly, Lz), psdSizes=psdSizes, psdCumm=psdCumm, periodic=True, seed=7)


# Insert the packing
sp.toSimulation(material="granr")

# Define engines
# more info: https://yade-dem.org/doc/introduction.html
O.engines = [
  # Reset forces
  ForceResetter(),
  # Approximate collision detection, create interactions
  InsertionSortCollider(
        # Sphere axis-aligned bounding box
        [Bo1_Sphere_Aabb()],
        label='collider'),
  # Handle interactions
  InteractionLoop(
        # Geom functors to create interaction geometry
        [Ig2_Sphere_Sphere_ScGeom()], 
        # Phys functors to create interaction from materials
        [Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
        # Law fcunctors to resolve interaction by computing forces per constituitive law
        [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'
  ),
  # Update positions uising Newton's equations
  NewtonIntegrator(
        # Non-viscous newton damping
        damping=.2,
        # Set gravity
        gravity=(0,0,0),
        # Create label for integrator engine
        label='newton'
  ),
  PeriTriaxController(
        # Create label/var name for triax engine
        label='triax',
        # Specify target values and whether they are strains or stresses
        goal=(sigma_iso, sigma_iso, sigma_iso),
        stressMask=7,
        # Type of servo-control
        dynCell=True,
        maxStrainRate=(max_strain_rate_iso, max_strain_rate_iso, max_strain_rate_iso),
        # Wait until the unbalanced force goes below this value
        maxUnbalanced=max_unbalanced_iso,
        relStressTol=rel_stress_tol_iso,
        # Call this function when goal is reached and the packing is stable
        doneHook='compactionFinished()'
  ),
  PyRunner(command='addPlotData()', iterPeriod=100),
]

# Add custom data to plot
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],
    n0=yade._utils.porosity(),
    #n1=yade._utils.voxelPorosity(200,*yade._utils.aabbExtrema()),
    n1=yade._utils.voxelPorosity(200,yade._utils.aabbExtrema()[0]*1.05,yade._utils.aabbExtrema()[0]*0.95),
    # mechanical coordination number Zm
    Zm=utils.avgNumInteractions(skipFree=True),
    # coordination number Cn
    Cn=utils.avgNumInteractions(skipFree=False),
    # 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'),
  '  i': ('Zm', None, 'n0', 'n1'),
  'ezz': ('sxx','syy','szz')
}
# Show the plot
plot.plot()


def compactionFinished():
  """
  Defines hook for setting up shearing stage after compaction is  finished.
  """
  # Reassign friction values for shearing stage
  # ..for future contacts change material
  O.materials[0].frictionAngle = radians(fric_final) # radians
  O.materials[0].etaRoll = eta_roll_final
  # ..for existing contacts, set contact friction directly
  for i in O.interactions:
    i.phys.tangensOfFrictionAngle = tan(radians(fric_final))
    i.phys.maxRollPl = min(eta_roll_final*O.bodies[i.id1].shape.radius, 
                           eta_roll_final*O.bodies[i.id2].shape.radius)  
  # 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 = (sigma_iso, sigma_iso, -strain_target)
  triax.stressMask = 3
  # Allow faster deformation along x,y to better maintain stresses
  triax.maxStrainRate = (max_strain_rate_horiz, max_strain_rate_horiz, max_strain_rate_vertical)
  # Next time, call triaxFinished instead of compactionFinished
  triax.doneHook = 'triaxFinished()'
  # Do not wait for stabilization before calling triaxFinished
  triax.maxUnbalanced = 10

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

O.run()

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