time for reaching a static equilibrium

Asked by Cloud on 2019-05-08

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:

from yade import pack
import math
nRead=utils.readParamsFromTable(
        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
)
from yade.params import table

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,
    frictionAngle=radians(compFricDegree),
    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
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

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

O.save('InitialState'+key+'.yade.gz')
print "### Initial state saved ###"

Thank you in advance!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Cloud
Solved:
2019-05-09
Last query:
2019-05-09
Last reply:
2019-05-08
Jan Stránský (honzik) said : #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 : #2

Hi Jan,

My bad.
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 : #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 : #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 : #5

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