No interactions while buiding small packing compared to the cell dimensions

Asked by Clerc Adriane

Hello,

I met a surprising behaviour while preparing a sample of spheres, could you tell me if there is a rational answer to it ?
With the script given in this message, I build a packing of 100 spheres in a cell of dimensions 0.15*0.15*0.15. Then in a second and third time I compact it with internal and external compaction.
During the 2nd phase, my sphere should swell and should enter in contact with each other and/or the walls. However, it seems that there is no interactions between spheres and between spheres and walls, the spheres just merge and their radius blow up...

The very interesting thing is that if I make the cell smaller (0.05*0.05*0.05 in line 21 of the script ) or if I increase the number of spheres (500 in line 14), then it works perfectly...

What could explain that ?

####SCRIPT####

from yade import pack, plot
from math import pi
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
load1=-20*10**3
load2=-100*10**3
young=356e6
poisson=0.42
frictAngle=32
density=2400
radiusRange=14.6/4.2
num_spheres=100
stabilityThreshold=0.015
EcTarget=10**-5
damp=0.05
# Create 6 bounding walls in order : left, right, bottom, top, back, front
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,
                            density=density,label='walls'))
dimCell=(0.15,0.15,0.15)

wall_ids=[]
wall_ids+=[O.bodies.append(utils.wall(0,axis=0,sense=1))] # left wall
wall_ids+=[O.bodies.append(utils.wall(dimCell[0],axis=0,sense=-1))] # right wall
wall_ids+=[O.bodies.append(utils.wall(0,axis=1,sense=1))] # bottom wall
wall_ids+=[O.bodies.append(utils.wall(dimCell[1],axis=1,sense=-1))] # top wall
wall_ids+=[O.bodies.append(utils.wall(0,axis=2,sense=1))] # back wall
wall_ids+=[O.bodies.append(utils.wall(dimCell[2],axis=2,sense=-1))] # front wall

# Create the sphere assembly
O.materials.append(FrictMat(young=young,poisson=poisson,
                            frictionAngle=radians(frictAngle),
                            density=density,label='spheres'))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),dimCell,0.3*(14.6+4.2)/4*10**-3,float(radiusRange-1)/(radiusRange+1),
             num_spheres,False,seed=0)
sp.toSimulation()

###############################
# Engine definitions
###############################

triax=TriaxialStressController(label='triax',
                height0=dimCell[1],width0=dimCell[0],depth0=dimCell[2],
                wall_left_id=wall_ids[0],wall_right_id=wall_ids[1],
                wall_bottom_id=wall_ids[2],wall_top_id=wall_ids[3],
                wall_back_id=wall_ids[4],wall_front_id=wall_ids[5],
                # Isotropic initial compression
                stressMask=7,
                goal1=load1,
                goal2=load1,
                goal3=load1,
                internalCompaction=True,
                # Comment if internalCompaction=True
                max_vel=0.01,stressDamping=0.25,strainDamping=0.8,
                # Uncomment if internalCompaction=True
                maxMultiplier=1.0001,
                finalMaxMultiplier=1.000001
                )

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
                 [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
                 [Ip2_FrictMat_FrictMat_FrictPhys()],
                 [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
# GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,
# timestepSafetyCoefficient=0.8),
        triax,
        NewtonIntegrator(damping=damp,label='newton')
        ]

#######################################################
# Start isotropic compression with internal compaction
#######################################################
O.dt=0.9*utils.PWaveTimeStep()

test=False
while test==False:
    O.run(1000,True)
    unb=unbalancedForce()
    mean=triax.meanStress
    test=unb<stabilityThreshold and abs((mean-load1)/load1)<0.01
    print '\r unbalanced force:',unb,' mean stress: ',triax.meanStress, O.iter,
    sys.stdout.flush()

O.save('confinedState.yade.gz')
print "\n ### Isotropic state saved ###"

#######################################################
# Start isotropic compression without internal compaction
#######################################################
triax.internalCompaction=False
triax.goal1=triax.goal2=triax.goal3=load2
# Enable energy tracking
O.trackEnergy=True

with open('record_Ec.txt','w') as f:
    f.write('Iter Ec\n')

test=False
while test==False:
    O.run(1000,True)
    Ec=O.energy['kinetic']
    mean=triax.meanStress
    test=Ec<EcTarget and abs((mean-load2)/load2)<0.01
    with open('record_Ec.txt','a') as f:
        f.write(str(O.iter)+' '+str(Ec)+'\n')
    print '\r Ec:',Ec,' mean stress: ',triax.meanStress, O.iter,
    sys.stdout.flush()

Question information

Language:
English Edit question
Status:
Expired
For:
Ubuntu yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Launchpad Janitor (janitor) said :
#1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.