yade-users team mailing list archive

Thread
Date
[Question #701581]: Difference between FrictMat and cohFrictMat without cohesion and moment

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Aoxi Zhang <question701581@xxxxxxxxxxxxxxxxxxxxx>
Date: Thu, 28 Apr 2022 14:05:51 -0000
Reply-to: question701581@xxxxxxxxxxxxxxxxxxxxx
Sender: noreply@xxxxxxxxxxxxx
New question #701581 on Yade:
https://answers.launchpad.net/yade/+question/701581

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



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