Facing problem with option in YADE that simulates grains similar to particle packing through Boundary contraction procedure in PFC 3D

Asked by Sarath Chandra Reddy on 2021-05-21

Dear friends,
I am a beginner in YADE. I am trying to simulate CD triaxial test. As a initial step, I am trying to create particle packing. I found that their are different ways for doing this i.e., using makeCloud, pack.randomDensePack with Compaction control and particle resizing.
I thick that pack.randomDensePack with Compaction control is similar to boundary contraction option in PFC 3D. But when I am trying to generate particle packing using this option, the unbalanced forces are not becoming less than the threshold value even when I tried different parameter values. Please advice on where I am going wrong with the parameter values or is there any mistakes in the code.

#Code for simple triaxial test
#Specify the variables
sigmaIso = 1.5e5 #Confining pressure
rMean = 0.00027 #Mean radius of the particles
rRelFuzz = 0.3 #The relative fuzz of the distribution
mn,mx = Vector3(0,0,0),Vector3(.00565,.00565,.00565)#The boundaries of the bounding box
initialFriction = 0.1 #Friction in radian
young=100e9 #Young's modulus of the particles
wallFriction=0 #Friction in the triaxial walls
num_spheres=6000#Predicted maximum number of spheres in simulation
rate = 1
targetPorosity=0.35 #Targeted porosity
finalFriction = .25 #Friction angle during deviatoric loading phase
sphere_density = 2650 #Density of the spheres
wall_density = 4000 #Density of the wall
wall_thickness = 0 #thickness of the walls

#Importing the relevant modules from yade
from yade import pack, qt, plot
O.materials.append(FrictMat(density=wall_density,young=young,frictionAngle=wallFriction, label='walls'))
O.materials.append(FrictMat(young=young,density=sphere_density,frictionAngle=initialFriction,label='sphere'))
#Appending the walls before we form the packing
walls=aabbWalls([mn,mx],thickness=wall_thickness, material='walls')
wallIds=O.bodies.append(walls)
#Using predicate to generate the sphere packing inside the walls
pred=pack.inAlignedBox((mn),(mx))
#Predicate and memorizeDb command saves it so that the packing can be used over and over again
sp=SpherePack()
sp=pack.randomDensePack(pred,radius=rMean,rRelFuzz=rRelFuzz,memoizeDb='/tmp/triaxpackcache.sqlite',returnSpherePack=True,material='spheres')
sp.toSimulation()
#Printing the mass of the spheres
print (utils.getSpheresMass(-1))
#Triaxial stress Controller
triax = TriaxialStressController(wall_bottom_id=wallIds[2],
                                 wall_top_id=wallIds[3],
                                 wall_left_id=wallIds[0],
                                 wall_right_id=wallIds[1],
                                 wall_back_id=wallIds[4],
                                 wall_front_id=wallIds[5],
                                 thickness=wall_thickness,
                                 #Stress controlled on three axes i.e., x*1+y*2+z*4
                                 stressMask=7,
                                 internalCompaction=False,
                                 label="triax"
                                 )
#Now we will start the engine
O.engines=[
           ForceResetter(),
           InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
           InteractionLoop(
           [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
           [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.1,es=0.1)], [Law2_ScGeom_MindlinPhys_Mindlin()]),
           GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
           triax,
           NewtonIntegrator(damping=0.2)
           ]
#Applying Confining Pressure
triax.goal1=triax.goal2=triax.goal3=-sigmaIso

while 1:
        O.run(1000,True)
        unb=unbalancedForce(),
        print ('unbalanced force:',unb,'mean stress:',triax.meanStress,'pressure:',abs(-sigmaIso-triax.meanStress)/sigmaIso)
        if unbalancedForce() < .01 and (abs(-sigmaIso-triax.meanStress)/sigmaIso) < .1:
               break

print ("Particle Distribution is Stable")

import sys

#while triax.porosity>targetPorosity:
 # initialFriction=0.95*initialFriction
 # setContactFriction(initialFriction)
  # print ("Friction:",initialFriction,"Porosity:",triax.porosity)
  # sys.stdout.flush()
   # O.run(500,1)
print ("Time for the deviator loads")

setContactFriction(finalFriction)
O.run(1000,True)
triax.internalCompaction=False
setContactFriction(finalFriction)
triax.stressmask=5
triax.goal2=-rate
triax.goal1=-sigmaIso
triax.goal2=-sigmaIso
O.engines = O.engines+[PyRunner(command='addPlotData()',iterPeriod=50000)]

def addPlotData():
# Lets start by defining a function that adds the data to a text file
    e11=-triax.strain[0]
    e22=-triax.strain[1]
    e33=-triax.strain[2]
    s11=-triax.stress(triax.wall_right_id)[0]
    s22=-triax.stress(triax.wall_top_id)[1]
    s33=-triax.stress(triax.wall_front_id)[2]
    plot.addData(Strain=e22,StressX=s11,StressY=s22,StressZ=s33,DeviatorStress=((s22-((s11+s33)/2))*1e-6))
    plot.saveDataTxt('data.txt',vars=('Strain','DeviatorStress','StressX','StressY','StressZ'))

plot.plots={'Strain':('DeviatorStress')}
plot.plot()

    #Finally we run the simulation automatically
O.run()

Thanking you in advance,
Sarath

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Sarath Chandra Reddy
Solved:
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#1

(sorry for sending earlier answer to users list, here it is)

Hi,
The problem is you parameters are extreme and the servo control of boundaries doesn't like that:
Super-high stiffness with super low confinement. If you visualize interaction forces in the 3D view you can see that they keep blinking, that's because the walls struggle to maintain a constant stress and the force oscillates between zero and something.

The particular case also leads automatic timestep to not be correct (may also be linked to the viscous Mindlin law you are using).
Without changing the physical problem I could stabilize your script with:
timestepSafetyCoefficient=0.25
triax.max_vel=0.002

I would discourage generating a randomDensePack for a triaxial loading since, what it does is actually another triaxial loading. Then when you load it as an initial state you have a discontinuity in mechanical properties which make it unstable again. Especially in your super-stiff case. Just a waste of computation time.
Better do everything in same conditions.

I would also avoid using super low confinement-to-stiffness ratio (1e-6 in your case) if possible, since it makes the whole thing not very efficient.

I'm surprized that randomDensePack seems to not reload previous packs despite you using memoizeDb. Do you see the same behaviour?

I hope it helps

Bruno

Revision history for this message
Sarath Chandra Reddy (sarath4321cityu) said :
#2

Thank you, Dr. Bruno, for the prompt reply and for sharing your knowledge of how the servo engine works.
The problem is your parameters are extreme and the servo control of boundaries doesn't like that:
Super-high stiffness with super low confinement.
Here, I used higher particle young’s modulus value as in the Mindlin model the stiffness also depends on contact overlap. As I understood previously from PFC 3D that even at E = 100 GPa, the normal stiffness values are in the range of 10^4 – 10^5 N/m after the confinement stage as compared to the Cundall and Stark model (Linear model) where people use normal stiffness in the range of 10^8 N/m.
If you visualize interaction forces in the 3D view you can see that they keep blinking, that's because the walls struggle to maintain constant stress and the force oscillates between zero and something.
Yes, I also observed this.
The particular case also leads automatic timestep to not be correct (may also be linked to the viscous Mindlin law you are using).
Without changing the physical problem I could stabilize your script with:
timestepSafetyCoefficient=0.25
triax.max_vel=0.002
I will try this
I would also avoid using a super low confinement-to-stiffness ratio (1e-6 in your case) if possible since it makes the whole thing not very efficient.
okay
I'm surprised that randomDensePack seems to not reload previous packs despite you using memoizeDb. Do you see the same behavior?
Yes, I was also surprized to see this behavior. This is one of the reasons, I raised the query as I doubted that I did something wrong with the code.

I would discourage generating a randomDensePack for a triaxial loading since, what it does is actually another triaxial loading. Then when you load it as an initial state you have a discontinuity in mechanical properties which makes it unstable again. Especially in your super-stiff case. Just a waste of computation time.
Better do everything in same conditions.

Are suggesting to use other ways for particle packing generation.
Can you please let me know if you recommend makeCloud with gravity settlement and using the Boundary contraction procedure be used to attain approximately the required porosity works better for this case? Please suggest if any other better way is available.

My study is more focused on the deviator stress vs strain behavior just after confinement. I wanted a mechanism that can replicate the effect of particle properties on initial slope of the deviator stress vs strain curve, i.e., with increase in particle young’s modulus, possion’s ratio, the slope of deviator stress versus strain should increase.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#3

Please check examples/triax-tutorial/script1.py, where the same engine is used to generate initial packing then apply deviatoric loading.
You can grow particle size or apply confinement by moving the walls.
Bruno

Revision history for this message
Sarath Chandra Reddy (sarath4321cityu) said :
#4

Thank you Dr.Bruno for the information. I will try this.