yade-users team mailing list archive

[Leonard <question692392@xxxxxxxxxxxxxxxxxxxxx>]

Question #692392 on Yade changed:
https://answers.launchpad.net/yade/+question/692392

    Status: Answered => Open

Leonard is still having a problem:
Hi Jan,
Many thanks for your reply.

With your suggestion, I got the i.geom.penetrationDepth=someValue. It is
really good compared to Nan before.

However, I found another two problems that I'd like to discuss:
(1) In the Script2(i.e., loading sample). No matter I set O.dt=0 or not, the i.geom.penetrationDepth (for a certain interaction) is the same. You may see it by skipping O.dt=0 in Script2 below.

(2) At the end of Script1 (i.e., making sample), I record a penetrationDepth of a certain interaction, for example O.interactions[2154,3930].geom.penetrationDepth=1.54203739497e-06. 
And when I go to Script2, I check the penetrationDepth of that interaction , the value I got in Script2 is slightly different with Script1 (i.e., O.interactions[2154,3930].geom.penetrationDepth=1.5419131404164482e-06). which means that the state has changed although it is slight. 

You may reproduce the problem using the following MWE (Script1 and
Script2) below if you are interested. (if O.interactions[2154,3930]
isn't existed, You may need to select a new one from script1):

##################Script 1#######################
from __future__ import division
from yade import pack, plot
num_spheres=7000#
targetPorosity = 0.45
compFricDegree = 30
finalFricDegree = 18
rate=-0.01
damp=0.6
stabilityThreshold=0.001
young=2e8
confinement=100e3

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

MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
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='spheres') for center,rad in sp])
triax=TriaxialStressController(
    maxMultiplier=1.+3e7/young,
    finalMaxMultiplier=1.+16e4/young,
    thickness = 0,
    stressMask = 7,
    internalCompaction=True,
)
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),
    triax,
    TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
    newton,
]

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

while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  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 "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

O.save('sample.yade.gz')
print "### Compacted state saved ###"

for i in O.interactions:
    print i.id1,i.id2,i.geom.penetrationDepth

##################Script 2#######################
from __future__ import division
from yade import pack, plot

O.load('sample.yade.gz')
confinement=100e3
triax=TriaxialStressController(
    thickness = 0,
    stressMask = 5,
    internalCompaction=False,
)
newton=NewtonIntegrator(damping=0.4)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
         Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'),
    newton,
]
triax.goal2=-1
triax.goal1=-confinement
triax.goal3=-confinement

print "penetrationDepth:", O.interactions[2154,3930].geom.penetrationDepth
print "O.dt:", O.dt
# O.dt=0 #(or skip this step for my question 1)
O.step()
print "penetrationDepth:", O.interactions[2154,3930].geom.penetrationDepth
#############################
Do you have any ideas about this?

Thanks
Leonard

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