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 unb0.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()
NewtonInt
