The spheres explodes at the first iter

it seems the modulus of the walls and particles are too different
 but when this is fixed, the problem still exists when O.run(3000,True)

And when the explode happens,

for ip in O.interactions:
 print O.bodies[ip.id1].state.pos
 print O.bodies[ip.id1].state.mass
 print O.bodies[ip.id2].state.pos
 print O.bodies[ip.id2].state.mass
 break

will print (a interaction between walls and sphere):
 Vector3(0,0,0.2)
0.0
Vector3(nan,nan,nan)
0.0

this problem is solved by change "O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,density=0,label='walls'))"
to

>>
"O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,density=2600,label='walls'))"

I have no idea....

If you do not assign any material to particles (like you did in
O.bodies.append(ymport.text('particles.txt'))
and/or
 sp=pack.SpherePack()
 sp.makeCloud(mn,mx,rmean2,rfuzz,numBall,False, 0.95,seed=1)
 sp.toSimulation()
), then the last O.materials.appended material is used for them. In your case with density 0 and therefore of course Nan after NewtonIntegrator. Specifying material should fix it:

O.bodies.append(ymport.text(... , material='spheres'))
sp.toSimulation(material='spheres')

cheers
Jan

Thanks Jan Stránský, that solved my question.

Hi Jan, I found the makeCloud function seems doesn't create non-overlapping cloud. The script I use is followed below:

#############################
from yade import ymport,pack,export
import sys
#surf=ymport.stl('disc.stl',fixed=False)
#O.bodies.append(surf)
key='a'
mv_x=200e-3
mv_y=200e-3
mv_z=200e-3
sxx=-2e5
targetPorosity = 0.50
rlow=3e-3
rhi=rlow*1.5
rmean=(rlow+rhi)/2.0
rfuzz=(rhi-rlow)/2.0/rmean
rmean2=rmean/2.0
numBall=int(mv_x*mv_y*mv_z/(4.0/3.0*pi*rmean**3)*(1-targetPorosity))
print numBall
O.materials.append(FrictMat(young=20e9,poisson=0.5,frictionAngle=0,density=2600,label='spheres'))
O.materials.append(FrictMat(young=20e9,poisson=0.5,frictionAngle=0,density=2600,label='walls'))

mn,mx=Vector3(-mv_x/2.0,-mv_y/2.0,mv_z),Vector3(mv_x/2.0,mv_y/2.0,0)
walls=aabbWalls([mn,mx],thickness=0,material="walls")
wallIds=O.bodies.append(walls)
center=Vector3(0,0,mv_z/2.0)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rmean2,rfuzz,numBall,False, 0.95,seed=1)
sp.toSimulation(material="spheres")
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.state.pos=(b.state.pos-center)*1+center
export.text('particles.txt',-1)

def getPoro():
 Vol=0.0
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   Vol+=o.shape.radius**3*4.0/3.0*pi
 return 1-Vol/mv_x/mv_y/mv_z
radiusMult = ( (1.0 -targetPorosity ) / (1.0 - getPoro()) )**(1.0/3.0)

print radiusMult

def expand(mult):
 print "expanding radius * %f" % mult
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   b.shape.radius*=mult
   b.state.mass=2600*pi*4.0/3.0*b.shape.radius**3
def calm():
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   b.state.vel=Vector3(0,0,0)
   b.state.angVel=Vector3(0,0,0)
print getPoro()
triax=ThreeDTriaxialEngine(
 maxMultiplier=1.+2.0e4/20e9,
 finalMaxMultiplier=1.+2.0e3/20e9,
 thickness = 0,
 stressControl_1 = True,
 stressControl_2 = True,
 stressControl_3 = True,
 ## Independant stress values for anisotropic loadings
 goal1=sxx,
 goal2=sxx,
 goal3=sxx,
 internalCompaction=True,
 Key=key,
)
newton=NewtonIntegrator(damping=0.7)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 newton,
 VTKRecorder(fileName='vtk/gravity-',recorders=['all'],label='vvtk',iterPeriod=1000)
]
O.dt=.5*PWaveTimeStep()
def equ(n):
 M=radiusMult**(1/float(n))
 for nn in range(0,n):
  expand(M)
  print 'step1: current poro: ', getPoro()
  for nn in range(0,50):
   O.run(5,True)
   calm()
  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
   calm
   if unb<1e-4:
    break
##################################################

#for ip in O.interactions:
# print O.bodies[ip.id1].state.pos
# print O.bodies[ip.id1].state.mass
# print O.bodies[ip.id2].state.pos
# print O.bodies[ip.id2].state.mass
# print ip.phys.kn
# print ip.phys.normalForce
# break
yade.qt.Controller(), yade.qt.View()
#equ(5)
O.run(1,True)
def hh():
 for ip in O.interactions:
  #print O.bodies[ip.id1].state.vel
  #print O.bodies[ip.id2].state.vel
  print ip.geom.penetrationDepth/rmean2
hh()
################################################################

the print output is usually 0.5-0.9, which is very large to make the particles explodes out of the box

>I found the makeCloud function seems doesn't create non-overlapping cloud

Not an option, no...
On the other hand I see that you implemented your own growing algorithm and it is incomplete.
I would suggest to check the function growParticles(), which also updates the interactions.
Bruno

Hi bruno,
 yes I'm using the way of expanding particles (usual way in PFC to generate dense particles with desired porosity). The first step is to create loose and non-overlapping particles. Unfortunately, sth seems wrong with makeCloud. Not sure if it is my own fault... When I export particles generated by PFC command (generate command, loose and non-overlapping particles before expanding their radius ), then the simulation goes well.

>sth seems wrong with makeCloud

Not an option, no...
On the other hand I see that you implemented your own growing algorithm and it is incomplete.
I would suggest to check the function growParticles(), which also updates the interactions.
Bruno

Hi bruno, the ThreeDTriaxialEngine is not activated in the above script, the problem may not be caused by the growParticles()....

I wonder how #7 and #9 can be so unclear. Let me try again.
1. makeCloud is not generating overlaps, definitely not. If you think it does please show a proof.
2. You growing method is wrong, you need to use growParticles() instead
(and better not ThreeDTriaxialEngine, TriaxialStressController is more flexible and more maintained)
Bruno

I incidentally found that my above statement was incorrect. "makeCloud is not generating overlaps" is only true under the assumption of consistent input by user.
It is admitted, for instance, that min should be smaller than max...
I hope it helps.
Bruno

Hi Bruno, I have solved the problems by fixing the min and max value, what a stupid mistake....

And besides the makeCloud function, the AABBWalls is also affected by this mistake.
Thank you for your patience