# time for reaching a static equilibrium

Hi,

I want to ask a question I have encountered, every time I change box size, it becomes very difficult to reach equilibrium. for example,mn,mx=Vector3(0,0,0),Vector3(1e-3,1e-3,1e-3) is Ok, but I change it to mn,mx=Vector3(0,0,0),Vector3(1e-2,1e-2,1e-2), it will take a long time or hard to reach equilibrium.

MWE:

import math
young = 0.858e9,
poisson = 0.3,
compFricDegree = 8,
finalFricDegree = 30.0,
num_spheres = 10000,
#young_wall=3.0e7
confiningS = 100e3,
isoconfining = 5e3,
AxialStrainLimit = 0.3,
rate = -0.02,
damp = 0.2,
unknownOk = True
)

mn,mx=Vector3(0,0,0),Vector3(1e-2,1e-2,1e-2) # corners of the initial packing
young = table.young
poisson = table.poisson
compFricDegree = table.compFricDegree
finalFricDegree = table.finalFricDegree
num_spheres = table.num_spheres
confiningS = table.confiningS
isoconfining = table.isoconfining
damp = table.damp
rate = table.rate
AxialStrainLimit=table.AxialStrainLimit
key='_cf-'+str(compFricDegree)+'_ff-'+str(finalFricDegree)+'_Y-'+str(young/1e6)+'_p-'+str(poisson)+'_S-'+str(confiningS/1000)
cd = 1 # cd=1 - consolidate drained cd=0- undrained
porosity = 0.72
stabilityThreshold = 0.001
thick = 0.01
d_limit = 0.15 # (d_limit/dmin = 0.25)
roll_stiff = 0.01
targetPorosity = 0.43

O.materials.append(FrictMat(
young=young,
poisson=poisson,
density=2650,
label='spheres')
)

O.materials.append(FrictMat(
young=young/10,
poisson=poisson,
frictionAngle=0,
density=0,
label='walls')
)
# create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,1.0/7.0,num_spheres,False,porosity,seed=1)# porosity=0.65 is the default value

triax=TriaxialStressController(
## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions. (Documention p290)
maxMultiplier=1.001, # spheres growing factor (fast growth)
finalMaxMultiplier=1.00001, # spheres growing factor (slow growth)
thickness = thick,
stressMask = 7, # 1*1 + 1*2 + 1*3 = 7 means x,y,z are on
internalCompaction=True, # if true the confining pressure is generated by growing particles.
)

newton = NewtonIntegrator(damping=damp)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys()],
[Law2_ScGeom_MindlinPhys_Mindlin()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
TriaxialStateRecorder(iterPeriod=1000,file='TriaxialRecorder'+key,truncate=1),
newton,
]

triax.goal1=triax.goal2=triax.goal3=-isoconfining #(5kpa)
timing.reset()
while 1:
O.run(1000, True)
unb=unbalancedForce()
print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
if unb < stabilityThreshold and abs(-isoconfining-triax.meanStress)/isoconfining<0.001:
break

print "### Initial state saved ###"

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Cloud
Solved:
2019-05-09
Last query:
2019-05-09
2019-05-08
 Jan Stránský (honzik) said on 2019-05-08: #1

Hello,

> I want to ask a question

what is the question (there is none in your post) ?

> mn,mx=Vector3(0,0,0),Vector3(1e-3,1e-3,1e-3) is Ok, but I change it to mn,mx=Vector3(0,0,0),Vector3(1e-2,1e-2,1e-2), it will take a long time

in the second case, the situation is 1000x larger (probably having 1000x more particles). No surprise it takes longer..
Also please be more specific what "long time" and "hard to reach" means. More real time? virtual time? it converges slower to equilibrium? does not converge at all? ...

cheers
Jan

 Cloud (strife-cloud) said on 2019-05-09: #2

Hi Jan,

My question is that Why is it very time-consuming to converge after the box size is modified to be larger ?(mn,mx=Vector3(0,0,0),Vector3(1e-3,1e-3,1e-3) -> Vector3(0,0,0),Vector3(1e-2,1e-2,1e-2)) .

>in the second case, the situation is 1000x larger (probably having 1000x more particles). No surprise it takes longer..

Yes, I agree. My setting is sp.makeCloud(mn,mx,-1,1.0/7.0,num_spheres=10000,False,porosity,seed=1),so the number of particles is 10000.

>Also please be more specific what "long time" and "hard to reach" means. More real time? virtual time? it converges slower to equilibrium? does not converge at all? ..

I have tried the following situations:

1) use GlobalStiffnessTimeStepper first: it's easy to converge to equilibrium when the box size is ((0,0,0),(1e-3,1e-3,1e-3)); it can not coverge at all when ((0,0,0),(1e-2,1e-2,1e-2)).(always unbalanced force: nan mean stress: nan)

2) use 0.5*utils.PWaveTimeStep(): it's also easy to converge to equilibrium when the box size is ((0,0,0),(1e-3,1e-3,1e-3)); it can not coverge at all when ((0,0,0),(1e-2,1e-2,1e-2)).(always unbalanced force: nan mean stress: nan)

3) try to set a timestep(for example O.dt = 1e-7 or some other values, i don't know how to set a proper value): something similar happened: smaller size is OK, bigger one takes a long real time to converge.

cheers,
Cloud

 Cloud (strife-cloud) said on 2019-05-09: #3

Beside, I found something interesting: I modify the box size in my MWE to Vector3(1e-3,1e-3,1e-3)
or Vector3(1,1,1) or Vector3(1e-1,1e-1,1e-1), all converge to equilibrium fast. It means only size of Vector3(1e-2,1e-2,1e-2) have problem of converge.

 Jan Stránský (honzik) said on 2019-05-09: #4

> always unbalanced force: nan mean stress: nan)

Good findings, now the right question is, why the unbalanced force is NaN. One situation is that there are no interactions. Also please try to find if some other value is NaN, like b.state.pos of some body or something like that.

There should not be much difference between PWaveTimeStep, GlobalStiffnessTimeStepper and setting O.dt directly (to a reasonable value)..

cheers
Jan

 Cloud (strife-cloud) said on 2019-05-09: #5

Hi, Jan, thank you for you response! I will try it.