non zero incident velocity on sticking contact

Hello,

please read [1] and provide more information and a MWE.

> I get ... even for non sliding
> ... get ... for sticking contacts

How exactly do you define "non sliding" and "sticking"?
How do you know that a contact is "even non sliding"?

> Why is that the case ...

Difficult to tell without more info (MWE).

> how can I get the 'correct' relative velocity at the contact

How do you define "correct relative velocity at the contact"?
How does it differ from the incidentVel approach?

> close to zero

How you define "close"?

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Hello again,

sorry, I did not think it would be helpful in this case, since this is about output. Here it is, modified from Bruno's triax example:

# -*- coding: utf-8 -*-
#*************************************************************************
# Copyright (C) 2010 by Bruno Chareyre *
# bruno.chareyre_at_grenoble-inp.fr *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************/

## This script details the simulation of a triaxial test on sphere packings using Yade
## See the associated pdf file for detailed exercises
## the algorithms presented here have been used in published papers, namely:
## * Chareyre et al. 2002 (http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
## * Chareyre and Villard 2005 (https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
## * Scholtès et al. 2009 (http://dx.doi.org/10.1016/j.ijengsci.2008.07.002)
## * Tong et al.2012 (http://dx.doi.org/10.2516/ogst/2012032)
##
## Most of the ideas were actually developped during my PhD.
## If you want to know more on micro-macro relations evaluated by triaxial simulations
## AND if you can read some french, it is here: http://tel.archives-ouvertes.fr/docs/00/48/68/07/PDF/Thesis.pdf

from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=500,# number of spheres
 compFricDegree = 30, # contact friction during the confining phase
 key='_triax_base_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.43 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
noub=0

def interaction2line(i):
    __x1,__x2,__x3=i.geom.contactPoint
    __nF1,__nF2,__nF3 = i.phys.normalForce
    __tF1,__tF2,__tF3 = i.phys.shearForce
    __v1,__v2,__v3= i.geom.incidentVel(i,avoidGranularRatcheting=True)
    fline = f"{i.id1:23} {i.id2:23} {__x1:23} {__x2:23} {__x3:23} {__v1:23} {__v2:23} {__v3:23} {__nF1:23} {__nF2:23} {__nF3:23} {__tF1:23} {__tF2:23} {__tF3:23}\n"
    return fline

def body2line(i):
    bline = f"\n"
    if isinstance (i.shape,Sphere):
        __idb = i.id
        __rad = i.shape.radius
        __p1,__p2,__p3=i.state.pos
        __v1,__v2,__v3=i.state.vel
        __w1,__w2,__w3=i.state.angVel
        bline = f"{__idb:23} {__rad:23} {__p1:23} {__p2:23} {__p3:23} {__v1:23} {__v2:23} {__v3:23} {__w1:23} {__w2:23} {__w3:23}\n"
    return bline

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE
## DEPENDING ON YOUR EDITOR, IT COULD BE DONE
## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES
## AND PRESSING CTRL+SHIFT+D

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break

##############################
### DEVIATORIC LOADING ###
##############################

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
#now goal2 is the target strain rate
triax.goal2=rate
# we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-10000
triax.goal3=-10000

#we can change damping here. What is the effect in your opinion?
newton.damping=0.1

#####################################################
### Example of how to record and plot data ###
#####################################################

#from yade import plot

def addPlotData():
    global noub
    namebod='b'+format(int(noub),'05d')+'.txt'
    namefor='f'+format(int(noub),'05d')+'.txt'
    lines = [body2line(i) for i in O.bodies]
    with open(namebod,"w") as f: # no need to close it if you use 'with'
            aline = f" iter {O.iter:23} time {O.time:23} porosity {utils.porosity():23}\n"
            f.writelines(aline)
            f.write(" body id body radius posx posy posz velx vely velz ang vel x ang vel y ang vel z '\n")
            f.writelines(lines)
    lines = [interaction2line(i) for i in O.interactions]
    with open(namefor,"w") as f: # no need to close it if you use 'with'
            f.write(" id1 id2 x1 x2 x3 v1 v2 v3 nF1 nF2 nF3 tF1 tF2 tF3 '\n")
            f.writelines(lines)
    noub+=1

# include a periodic engine calling that function in the simulation loop
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='addPlotData()',label='recorder')]+O.engines[5:7]

O.run(100,True)

'non sliding' or 'sticking' contacts (not a proper term, I know) are for me the ones where the ratio of the norm of the tangent to the normal force is much smaller than the tangent of the friction radius. I would expect no sliding at such a contact.

The relative velocity I get from the body velocities is exactly the same as the one returned by incidentVel. For contacts where I would expect no relative velocity, it is of the order of magnitude of the velocities of the bodies, here 10^{-3}. I would expect it to be significantly smaller, even if not zero.

What am I getting wrong?

Cheers,
eleni

> I did not think it would be helpful in this case, since this is about output. Here it is, modified from Bruno's triax example:

such code is indeed not much helpful :-(
What I asked is a MWE [1], M=minimal.
For this case it could be (by words):
I have these two spheres, I set this and that velocities and run one step. Velocities are this and that, forces are this and that, incidentVel is this but I expect that or zero.

> 'non sliding' or 'sticking' contacts ... are for me the ones ...

a "naive" definition could be what? Fs<<<Fn? Fn<<<Fs?

> the tangent of the friction radius

tangent of friction angle?

> The relative velocity I get from the body velocities is exactly the same as the one returned by incidentVel.

Thanks for clarification.

> I would expect no sliding at such a contact.
> For contacts where I would expect no relative velocity, it is of the order of magnitude of the velocities of the bodies, here 10^{-3}. I would expect it to be significantly smaller, even if not zero.

Why do you expect this?

In Yade, there are a few cases where you can force no sliding:
- clumps (rigid body set of a few particles)
- fixing ("blocking DOFs") of certain particles and prescribe their motion such that they move like rigid body

Otherwise, it is not possible to force two particles to "stick" or "no sliding".
And actually it is even not reasonable.
Image two "sticking" particles A and B. Now a particle C impacts to B, which would normally make it rotate.
Now should B rotate due to the impact? Should it do nothing doe to the sticking with A? ...

One of the basic concepts of Yade is the "simulation loop" [2]. In this case there are important the blue and the green part.
1) Based on position and velocities, forces are computed.
2) Based on forces and current velocity and position, new velocity and position is computed
1) ...
2) ...
...

So one interaction only partially influences its bodies motion. Body motion is a result of many interactions.
Maybe this may be source of the discrepancy of expectation and results?

Also, the "sticking" contacts according to your "ratio definition" is no special from the loop point of view, just the forces has some different values than in other cases. No "simulation loop" reason for the contact to be non sliding..

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/introduction.html#function-components

Hi again,

and thanks for taking the time.

As a minimum example: two spheres from the above simulation are in contact. The ratio of tangent to normal force is 0.1983, the friction coefficient is 0.5774. I would expect incidentVel for the contact to give almost zero, but it gives 6.8847e-04, one order of magnitude less than the linear velocities of the spheres in contact. Is it normal that there is relative shear displacement at a contact that does not satisfy the sliding criterion?

Cheers,
eleni

A clarification: I am not trying to force spheres not to slide, I just want to see if they do. According to the contact force they don't slide, according to the relative velocity at the contact they do. I don't understand why.

Hi,

particles have some contact stiffness both in normal and tangential direction. So they can have relative shear displacement at contact. The shear force is computed based on this displacement.

Cheers,
Karol

> I would expect incidentVel for the contact to give almost zero

why almost zero, why not strict zero? :-)
If almost zero, what is the limit where it is almost and what not?

> I would expect incidentVel for the contact to give almost zero, but it gives 6.8847e-04, one order of magnitude less than the linear velocities of the spheres in contact.

incidentVel is independent of absolute values of velocities of the particles. Only their difference matters.
You could have two particles at the speed of light (or even more in a model :-), still they can have arbitrary incidentVel (depending also significantly on rotation).
You could have two particles with fixed position, still they can have arbitrary incidentVel (bacause of rotation).

> Is it normal that there is relative shear displacement at a contact that does not satisfy the sliding criterion?

yes, absolutely

> According to the contact force they don't slide, according to the relative velocity at the contact they do.

depends on the definition of "slide".
They do not slide from friction point of view (there is no energy dissipation by friction).
But they still have mutual shear displacement/velocity, although in an elastic regime Fs=G*us, as explained by Karol.

As discussed, (under general loading conditions, not like nothing happens or free fall etc.) there is no reason why two "free" (non-clump etc.) spheres should have no contact "sliding" (zero mutual displacement/velocity, normal or shear).

The contact geometry (sliding, penetrationDepth etc.) influences (through computed forces) PARTIALLY the motion of particles. Actual motion of particles is SUM of ALL contacts.
The motion of particles FULLY determines the contact geometry (sliding etc.).
Therefore, you cannot expect certain exact value of contact geometry based only on that one contact.

Cheers
Jan

Thanks Karol, that answers it!