Difference between FrictMat and cohFrictMat without cohesion and moment

Asked by Aoxi Zhang

Hi,
I'd like to ask that is it the same between FrictMat and cohFrictMat without moment and cohesion?

As described in Yade documentation, cohFrictMat is the material extending FrictMat with cohesive properties and rotational stiffnesses, I am curious that if I close cohesion and moment in cohFrictMat, and use the same parameters as that of FrictMat, is the packing expected to be the same? (The corresponding O.engines-InteractionLoop should be modified as well)

The background of this question is: I got a problem that particle velocity varies with time step (asked at thread[1]), in which the material is cohFrictMat. Then I made an attempt to use FrictMat and went through exactly the same procedure (described below), and I found that the problem didn't occure in FrictMat case. Thereby, I guess the problem (particle velocity varies with time step) might be related to the material/ contact law I uesd.

The MWE for reproducing the question is attached below, it contains 3 phases (I know it is better to integrate them in 1 script, I am sorry that I didn't successfully make them in one go) :

(1) phase1_frict.py/ phase1_cohFrict.py is for sample generation (with frictMat/ cohFrictMat but same material parameters) and stress relaxation, it takes around 12 minutes if you want to reproduce the results.

The idea is: firstly having a sample according to Bruno's triaxial example, then closing triax engine and increasing damping to quickly reduce unb. phase1 stops when unb reaches a small value (e.g.,1e-7).

(2) phase2.py is the test under different dt, if you want to reproduce the resutls, you may use yade-batch mode, it takes around 11 minutes in general.

The idea of this test is: firstly load the sample generated from phase1, then turn damping to 0, close time stepper and set O.dt with different values, select some particles whose velocity will be recorded. phase2 will terminate when the time reaches 3e-4.

(3) phase3.py is for plotting the results to get a figure which compared the results under different dt.
##################### MWE:
#### phase1_frict.py######
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
from random import gauss
import timeit

num_spheres=5000
targetPorosity = 0.405
compFricDegree = 30
finalFricDegree = 19
rate=-0.01
damp=0.6
stabilityThreshold=0.001

confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=2e8,poisson=0.3,frictionAngle=0,density=0,label='walls'))

MatSand = O.materials.append(FrictMat(young=2e8,poisson=0.3,frictionAngle=radians(30),\
           density=2650.0,label='sand'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,label="TimeStepper"),
 TriaxialStressController(
  maxMultiplier=1.0+3e7/2e8,
  finalMaxMultiplier=1.+16e4/2e8,
  thickness = 0,
  stressMask = 7,
  internalCompaction=True,
  label='triax'),
 newton,
]

Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
   break

print "### state 1 Isotropic completed ###"
triax.internalCompaction=False
import sys
while triax.porosity - targetPorosity>0.0001:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

print "### state 2 Reach target porosity completed ###"
setContactFriction(radians(finalFricDegree))

NewtonIntegrator.damping=0.9
newton.damping=0.9 ## just make sure damping has been changed, according to Yade question:https://answers.launchpad.net/yade/+question/701438
O.engines=O.engines+[PyRunner(command='stopStressRelaxation()', iterPeriod=1,label="StopRelaxation")]

for i in range(6):
 O.bodies[i].state.vel = Vector3(0,0,0)
 O.bodies[i].state.blockedDOFs='xyzXYZ'

triax.dead=True

def stopStressRelaxation():
 print('unb',unbalancedForce())
 if unbalancedForce()<1e-7:
  O.pause()
  NewtonIntegrator.damping=0.0#newton.damping=0
  print "stress relaxation finished, damping has been set as 0"
  print "Please go to phase2"
  O.save('sample_generatedFromPhase1_afterRelax.yade.gz')

O.run()
utils.waitIfBatch()
############################
############ phase1_cohFrict.py###########
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
from random import gauss
import timeit
import pickle

num_spheres=5000
targetPorosity = 0.405
compFricDegree = 30
finalFricDegree = 19
rate=-0.01
damp=0.6
stabilityThreshold=0.001

confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=2e8,poisson=0.3,frictionAngle=0,density=0,label='walls'))

MatSand = O.materials.append(CohFrictMat(isCohesive=False,young=2e8,poisson=0.3,frictionAngle=radians(30),\
           density=2650.0,normalCohesion=0, shearCohesion=0,\
           momentRotationLaw=False,label='sand'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="Ip2Coh"),
   Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,label="TimeStepper"),
 TriaxialStressController(
  maxMultiplier=1.2,
  finalMaxMultiplier=1.004,
  thickness = 0,
  stressMask = 7,
  internalCompaction=True,
  label='triax'),
 newton,
]

Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
   break

print "### state 1 Isotropic completed ###"
triax.internalCompaction=False
import sys
while triax.porosity - targetPorosity>0.0001:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

print "### state 2 Reach target porosity completed ###"
setContactFriction(radians(finalFricDegree))

NewtonIntegrator.damping=0.9
newton.damping=0.9 ## just make sure damping has been changed, or see Yade question:https://answers.launchpad.net/yade/+question/701438
O.engines=O.engines+[PyRunner(command='stopStressRelaxation()', iterPeriod=1,label="StopRelaxation")]

for i in range(6):
 O.bodies[i].state.vel = Vector3(0,0,0)
 O.bodies[i].state.blockedDOFs='xyzXYZ'

triax.dead=True

def stopStressRelaxation():
 print('unb',unbalancedForce())
 if unbalancedForce()<1e-7:
  O.pause()
  NewtonIntegrator.damping=0.0#newton.damping=0
  print "stress relaxation finished, damping has been set as 0"
  print "Please go to phase2"
  O.save('sample_generatedFromPhase1_afterRelax.yade.gz')

O.run()
utils.waitIfBatch()
############################
########## phase2.py######
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import random
import pickle
from yade import ymport,export

utils.readParamsFromTable(dt=1e-7)
from yade.params import table

finalFricDegree=19

O.load('sample_generatedFromPhase1_afterRelax.yade.gz')

RunTimeLength=3e-4 # how long for the time of data recording
Odt=table.dt
O.dt=Odt

O.engines[3].dead=True## turn off GlobalStiffnessTimeStepper
# TimeStepper.dead=True## Note, use label cannot successfully close time stepper
StopRelaxation.dead=True## turn off StopRelaxation in previous phase
NewtonIntegrator.damping=0.0
newton.damping=0.0 ## just make sure damping has been changed, or see Yade question:https://answers.launchpad.net/yade/+question/701438
monitoredSand=[]
def getMonitoredSand():
 halfHeigh=aabbDim()[1]/2
 for i in O.bodies:
  if isinstance(i.shape, Sphere):
   if halfHeigh<i.state.pos[1]<halfHeigh+0.0015:
    i.shape.color=[1,0,0]
    monitoredSand.append(i.id)
   else:
    i.shape.color=[1,1,1]
 print ("number of monitoredSand is",len(monitoredSand))

getMonitoredSand()

def avgVel(idList):
 vel=0.0
 avg=0.0
 for i in idList:
   vel+=O.bodies[i].state.vel[1]
 avg=vel/len(idList)
 return avg

print("O.dt=",O.dt)
O.engines=O.engines+[PyRunner(command='outputData()', virtPeriod=RunTimeLength/1000,label="outputData")]
O.engines=O.engines+[PyRunner(command='stopSimulation()', iterPeriod=10,label="stopSimu")]

def stopSimulation():
 print ('Running, O.time,',O.time)
 if O.time > RunTimeLength:
  O.pause()
  print ('Running finished, O.realtime,',O.realtime)
  # plot.saveDataTxt('historyData_calmdownVel{}.txt'.format(calmDownVel)) # should revise the monitored data at history

def outputData():
 f = open("ResultsOfPhase2_dt{}.txt".format(Odt), "a+")
 f.write('{} {} {}\n'.format(O.time,O.dt,avgVel(monitoredSand)))
 f.close()

O.resetTime()
O.run()
utils.waitIfBatch()
#####################

###### phase3.py#####
from yade import pack, plot
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
import matplotlib.ticker as ticker
import numpy as np
import pandas as pd
font1 = {'family': 'Arial', 'fontweight': 'normal', 'size': 14, 'style':'normal'}

linewidth=2
fig , ax = plt.subplots(figsize=(6.25,3))
bwith = 0.5# here control the width of axis
ax.spines['bottom'].set_linewidth(bwith)
# ax.spines['bottom'].set_color('black')
ax.spines['left'].set_linewidth(bwith)
# ax.spines['left'].set_color('black')
ax.spines['top'].set_linewidth(bwith)
# ax.spines['top'].set_color('black')
ax.spines['right'].set_linewidth(bwith)
# ax.spines['right'].set_color('black')

def plotResults(dt):
 f=pd.read_csv("ResultsOfPhase2_dt{}.txt".format(dt),sep=" ",header=None)
 f.columns=["Otime","Odt","avgVel"]
 Otime,Odt,avgVel= np.array(f["Otime"]),np.array(f["Odt"]),np.array(f["avgVel"])
 ax.plot(Otime,avgVel,label="dt={}".format(dt))

plotResults(1e-7)
plotResults(1e-8)
plotResults(1e-9)
plotResults(5e-9)
# =============================================================================
# Figure setting
#1 set legend
ax.legend(loc=0,frameon=False,prop={'family': 'Arial','size': 11})

ax.set_xlabel('Time (s)',font1)
ax.set_ylabel('Velocity (m/s)',font1)

ax.tick_params(which='both',labelsize=10,direction='in')
plt.xticks(fontproperties = 'Arial', size = 10)
plt.yticks(fontproperties = 'Arial', size = 10)
# plt.ylim(-3e-5,-2e-5)
plt.xlim(0,3e-4)

plt.tight_layout(pad=0.4)
plt.savefig('Fig_compare_dt.jpg', dpi=600)
plt.show()
#####

#### For running phase2.py in batch mode, another txt file could be as below:
dt
1e-7
5e-8
1e-8
5e-9
1e-9
######################################################################

For the purpose of saving your time, if you would like to directly see the results, I also put all of them in [2].

Do you have any idea of why the problem (particle velocity varies with time step) only occurs in the case of cohFrictMat? In addition, why cohFrictMat without cohesion and moment doesn't behave the same as the case of FrictMat even with same input parameters?

Thanks!
Aoxi

[1]https://answers.launchpad.net/yade/+question/701438
[2]https://we.tl/t-Sam5oT0Ail

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

I think we can find logical that CohFrictMat without cohesion and moment would indeed be equivalent to FrictMat. Whether that is the case may indeed deserve attention, since the C++ implementations [1] look to be completely independent (no one is calling the other one, for instance).

If you seek help in investigating this question, I would advice you further try to propose minimal scripts, though.

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/CohesiveFrictionalContactLaw.cpp#L109 and https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ElasticContactLaw.cpp#L56

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#2

Hi,
Unless you found a bug, yes, the two cases should be equivalent.
It's hard to know if you made a mistake in converting materials or if there is really a difference due to the code.

Currently, your scripts are way too complex to dive in. If you could elaborate a simpler example maybe someone could check in more details.

Bruno

Can you help with this problem?

Provide an answer of your own, or ask Aoxi Zhang for more information if necessary.

To post a message you must log in.