The uniaxial compression experiment fails to reach the critical state

Asked by Zhicheng Gao

After isotropic consolidation, the confining pressure remains unchanged in the X and Y directions, and the axial compression in the Z-axis direction is 40%, which has not yet reached the critical state. Usually, the critical state can be reached at 25% of the axial strain. Why is this? The procedure is as follows:

from __future__ import print_function
from yade import pack, qt, plot
from math import *

nRead=readParamsFromTable(
        ## model parameters
        num_spheres=30000,
        targetPorosity= .32,
        confiningPressure=-50000,
        ## material parameters
        compFricDegree=25,#contact friction during the confining phase
        finalFricDegree=30,#contact friction during the deviatoric loading
        young=2e8,
        poisson=.2,
        density=2650,
        alphaKr=7.5,
        alphaKtw=0,
 competaRoll=.1,
        finaletaRoll=.1,
        etaTwist=0,
        normalCohesion=0,
        shearCohesion=0,
        ## fluid parameters
        fluidDensity=1000,
        dynamicViscosity=.001,
        ## control parameters
        damp=0,
        stabilityThreshold=.001,
        ## output specifications
        filename='suffusion',
        unknowOk=True
)

from yade.params.table import *

O.load('compactedStatesuffusion25.yade.gz')
##__________________________________third section, deviatoric loading__________________________

triax.width0=aabbDim()[0]
triax.height0=aabbDim()[1]
triax.depth0=aabbDim()[2]

# change control type: keep constant confinement in x,y, 40% compression in z
triax.stressMask=3
triax.goal3=-.1
triax.goal1=confiningPressure
triax.goal2=confiningPressure

# enable energy tracking in the code
O.trackEnergy=True
# define function to record history
def history():
    sig_xx=-triax.stress(triax.wall_right_id)[0]
    sig_yy=-triax.stress(triax.wall_top_id)[1]
    sig_zz=-triax.stress(triax.wall_front_id)[2]
    plot.addData(unbalanced=unbalancedForce(),i=O.iter,exx=-triax.strain[0],
            eyy=-triax.strain[1], ezz=-triax.strain[2],
            sxx=sig_xx,syy=sig_yy,szz=sig_zz,
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            dsratio=(-sig_zz+sig_xx)/((-sig_xx-sig_yy-sig_zz)/3.0),
            porosity=porosity(),Etot=O.energy.total(),**O.energy# save all available energy data
            )
O.engines=O.engines+[PyRunner(command='history()',iterPeriod=500)]
t0=O.time
while True:
    O.run(1000,True)
    print('running',-triax.strain[2])
    unb=unbalancedForce()
    if unb<stabilityThreshold and abs(-.4-triax.strain[2])/0.4<0.001:
        break
# define what to plot
plot.plots={'i':('unbalanced','porosity'),' i':('sxx','syy','szz',None,'exx','eyy','ezz'),' ezz':('dsratio'),'ezz':('ev')}
# show the plot
plot.plot()
# save plot data
plot.saveDataTxt('uneroded_soil.txt')

The calculation result is shown in the figure in the link:
 https://ibb.co/Yc7jbbQ

Obviously the volumetric strain is still changing and has not reached the critical state

Question information

Language:
English Edit question
Status:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

> The uniaxial compression ...
> After isotropic consolidation, the confining pressure remains unchanged in the X and Y directions ...

usually "uniaxial compression", uniaxial stress etc. refers to a scenario, where the lateral stress is zero.
This I would call triaxial compression.

> Usually, the critical state can be reached at 25% of the axial strain.

Please be more specific.
What is "usually"? In experiments? simulations?

> Why is this?

Because of the combination of packing, material model and loading :-)

Cheers
Jan

Revision history for this message
Robert Caulk (rcaulk) said :
#2

Hello,

Please, no external links [1].

>>Usually, the critical state can be reached at 25% of the axial strain. Why is this?

What is "usually?" Does that mean you are comparing to a common "triaxial" lab test with flexible membrane? Keep in mind that the walls are stiff and frictional in your simulation here so you should expect different behavior.

Cheers,

Robert

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#3

Hi, Jan
> Usually, the critical state can be reached at 25% of the axial strain.

This conclusion was obtained in many experiments and simulations。

> Because of the combination of packing, material model and loading :-)

Is this result normal?

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#4

Hi, Robert

>What is "usually?" Does that mean you are comparing to a common "triaxial" lab test with flexible membrane? Keep in mind that the walls are stiff and frictional in your simulation here so you should expect different behavior.

"Usually" means that many experiments and simulations have reached this conclusion, and they also use rigid walls. In my simulation the wall is frictionless.

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#5

Is it because my loading speed is too fast, but at present, it takes two days to calculate 40% axial strain, and the CPU is R7-5800H. Is there a problem with my code? PFC calculation is even faster than this

Revision history for this message
Robert Caulk (rcaulk) said :
#6

Hello,

Can you please provide a citation to an experiment with the same loading
conditions as you have here, plus the frictionless rigid walls that you are
using?

Cheers,

Robert

Le sam. 6 nov. 2021 à 15:11, Zhicheng Gao <
<email address hidden>> a écrit :

> Question #699362 on Yade changed:
> https://answers.launchpad.net/yade/+question/699362
>
> Zhicheng Gao gave more information on the question:
> Hi, Robert
>
> >What is "usually?" Does that mean you are comparing to a common
> "triaxial" lab test with flexible membrane? Keep in mind that the walls
> are stiff and frictional in your simulation here so you should expect
> different behavior.
>
> "Usually" means that many experiments and simulations have reached this
> conclusion, and they also use rigid walls. In my simulation the wall is
> frictionless.
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Jérôme Duriez (jduriez) said :
#7

Hello,

> PFC calculation is even faster than this

Then, why did you choose to use YADE ? ;-)

Revision history for this message
Jan Stránský (honzik) said :
#8

> Is there a problem with my code?

please provide the complete code if you want a serious answer

Cheers
Jan

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#9

Dear Jan,
 the code of isotropic compression is as follows:
##______________ First section, generate sample_________

from __future__ import print_function
from yade import pack, qt, plot
from math import *
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

nRead=readParamsFromTable(
        ## model parameters
        num_spheres=50000,
        targetPorosity= .32,
        confiningPressure=-50000,
        ## material parameters
        compFricDegree=30,#contact friction during the confining phase
        finalFricDegree=30,#contact friction during the deviatoric loading
        young=2e8,
        poisson=.2,
        density=2650,
        alphaKr=7.5,
        alphaKtw=0,
 competaRoll=.1,
        finaletaRoll=.1,
        etaTwist=0,
        normalCohesion=0,
        shearCohesion=0,
        ## fluid parameters
        fluidDensity=1000,
        dynamicViscosity=.001,
        ## control parameters
        damp=0,
        stabilityThreshold=.001,
        ## output specifications
        filename='suffusion',
        unknowOk=True
)

from yade.params.table import *

mn,mx=Vector3(0,0,0),Vector3(0.04,0.04,0.04)
psdSizes=[0.00042,0.0005,0.00208,0.0024]
psdCumm=[0.0,0.25,0.25,1.0]
# create materials for spheres
#shear strength is the sum of friction and adhesion, so the momentRotationLaw=True
O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
# generate particles packing
sp=pack.SpherePack()
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=1)
sp.toSimulation(material='spheres')

totalnumber=0
for b in O.bodies:
    if (isinstance(b.shape,Sphere)):
        totalnumber=totalnumber+1
print(totalnumber)

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)],
            [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)],
 ),
        #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        TriaxialStressController(label='triax',
            # specify target values and whether they are strains or stresses
            goal1=confiningPressure,goal2=confiningPressure,goal3=confiningPressure, stressMask=7,
            # type of servo-control, external walls compaction
            internalCompaction=False,
            thickness=0,
            ),
        NewtonIntegrator(damping=0.3)
]

qt.View()

O.dt=0.5*PWaveTimeStep()
import sys

def history():
    plot.addData(unbalanced=unbalancedForce(),i=O.iter,exx=-triax.strain[0],
            eyy=-triax.strain[1], ezz=-triax.strain[2],
            sxx=-triax.stress(triax.wall_right_id)[0],syy=-triax.stress(triax.wall_top_id)[1],szz=-triax.stress(triax.wall_front_id)[2],
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            porosity=porosity(),Etot=O.energy.total(),**O.energy# save all available energy data
            )
O.engines=O.engines+[PyRunner(command='history()',iterPeriod=200)]

while True:
    O.run(1000,True)
    unb=unbalancedForce()
    sxx=triax.stress(triax.wall_left_id)[0]
    syy=triax.stress(triax.wall_bottom_id)[1]
    szz=triax.stress(triax.wall_back_id)[2]
    print('unbalanced force:',unb,sxx,syy,szz)
    if unb<stabilityThreshold and abs((confiningPressure-sxx)/confiningPressure)<0.0001 and\
                                  abs((confiningPressure-syy)/confiningPressure)<0.0001 and\
                                  abs((confiningPressure-szz)/confiningPressure)<0.0001:
        break

print('Friction:',compFricDegree,'porosity:',triax.porosity)
# reaching a specified porosity precisely
while triax.porosity>targetPorosity:
    compFricDegree=0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print('Friction:',compFricDegree,'porosity:',triax.porosity)
    sys.stdout.flush()
    while True:
        O.run(500,True)
        unb=unbalancedForce()
        if unb<stabilityThreshold and abs((confiningPressure-sxx)/confiningPressure)<0.0001 and\
                                      abs((confiningPressure-syy)/confiningPressure)<0.0001 and\
                                      abs((confiningPressure-szz)/confiningPressure)<0.0001:
            break

# check the position of spheres
maxX=O.bodies[1].state.pos[0]
minX=O.bodies[0].state.pos[0]
maxY=O.bodies[3].state.pos[1]
minY=O.bodies[2].state.pos[1]
maxZ=O.bodies[5].state.pos[2]
minZ=O.bodies[4].state.pos[2]
ball_out_of_wall=list()
totalnumber2=0
for b in O.bodies:
    if (isinstance(b.shape,Sphere)):
        totalnumber2=totalnumber2+1
    if b.state.pos[0]>maxX or b.state.pos[0]<minX:
        ball_out_of_wall.append(b.id)
        print('the walls or balls are moving too fast, the time step is too big')
    elif b.state.pos[1]>maxY or b.state.pos[1]<minY:
        ball_out_of_wall.append(b.id)
        print('the walls or balls are moving too fast, the time step is too big')
    if b.state.pos[2]>maxZ or b.state.pos[2]<minZ:
        ball_out_of_wall.append(b.id)
        print('the walls or balls are moving too fast, the time step is too big')
print(ball_out_of_wall,totalnumber2)

# change the contact parameters to the final calibration value
setContactFriction(radians(finalFricDegree))
for b in O.bodies:
    b.mat.etaRoll=finaletaRoll
for i in O.interactions:
    i.phys.etaRoll=finaletaRoll

O.run(100000,1)
while True:
    O.run(1000,True)
    unb=unbalancedForce()
    sxx=triax.stress(triax.wall_left_id)[0]
    syy=triax.stress(triax.wall_bottom_id)[1]
    szz=triax.stress(triax.wall_back_id)[2]
    print('unbalanced force:',unb,sxx,syy,szz)
    if unb<stabilityThreshold and abs((confiningPressure-sxx)/confiningPressure)<0.0001 and\
                                  abs((confiningPressure-syy)/confiningPressure)<0.0001 and\
                                  abs((confiningPressure-szz)/confiningPressure)<0.0001:
        break
print(porosity())
O.save('compactedState'+filename+'25.yade.gz')
print('Compacted state saved')
binsSizes, binsProc, binsSumCum=psd(bins=20,mass=True)
pylab.semilogx(binsSizes, binsProc,label='Mass PSD of (free) %d random spheres'%len(binsSizes))
#pylab.show()
pylab.savefig('con.jpg')

# define what to plot
plot.plots={'i':('unbalanced','porosity'),' i':('sxx','syy','szz',None,'exx','eyy','ezz'),'ezz':('ev')}
# show the plot
plot.plot()

the code of triaxial compression is as follows:
from __future__ import print_function
from yade import pack, qt, plot
from math import *

nRead=readParamsFromTable(
        ## model parameters
        num_spheres=50000,
        targetPorosity= .32,
        confiningPressure=-50000,
        ## material parameters
        compFricDegree=0,#contact friction during the confining phase
        finalFricDegree=30,#contact friction during the deviatoric loading
        young=2e8,
        poisson=.2,
        density=2650,
        alphaKr=7.5,
        alphaKtw=0,
 competaRoll=.1,
        finaletaRoll=.1,
        etaTwist=0,
        normalCohesion=0,
        shearCohesion=0,
        ## fluid parameters
        fluidDensity=1000,
        dynamicViscosity=.001,
        ## control parameters
        damp=0,
        stabilityThreshold=.001,
        ## output specifications
        filename='suffusion',
        unknowOk=True
)

from yade.params.table import *

O.load('compactedStatesuffusion25.yade.gz')
##__________________________________third section, deviatoric loading__________________________
#triax.computeStressStrainInterval=1
triax.width0=aabbDim()[0]
triax.height0=aabbDim()[1]
triax.depth0=aabbDim()[2]
print(-triax.strain[0],-triax.strain[1],-triax.strain[2],O.iter)

# change control type: keep constant confinement in x,y, 40% compression in z
triax.stressMask=3
triax.goal3=-.1
triax.goal1=confiningPressure
triax.goal2=confiningPressure

step=O.iter
# enable energy tracking in the code
O.trackEnergy=True
# define function to record history
def history():
    sig_xx=-triax.stress(triax.wall_right_id)[0]
    sig_yy=-triax.stress(triax.wall_top_id)[1]
    sig_zz=-triax.stress(triax.wall_front_id)[2]
    ii=O.iter-step
    plot.addData(unbalanced=unbalancedForce(),i=ii,exx=-triax.strain[0],
            eyy=-triax.strain[1], ezz=-triax.strain[2],
            sxx=sig_xx,syy=sig_yy,szz=sig_zz,
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            dsratio=(-sig_zz+sig_xx)/((-sig_xx-sig_yy-sig_zz)/3.0),
            porosity=porosity(),Etot=O.energy.total(),**O.energy# save all available energy data
            )
O.engines=O.engines+[PyRunner(command='history()',iterPeriod=500)]

while True:
    O.run(1000,True)
    print('running',-triax.strain[2])
    unb=unbalancedForce()
    if unb<stabilityThreshold and abs(-.25-triax.strain[2])/.25<0.001:
        break
# define what to plot
plot.plots={'i':('unbalanced','porosity'),' i':('sxx','syy','szz',None,'exx','eyy','ezz'),' ezz':('dsratio'),'ezz':('ev')}
# show the plot
plot.plot()
# save plot data
plot.saveDataTxt('uneroded_soil.txt')
O.save('devia'+filename+'25.yade.gz')

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#10

Dear Jerome,
I used Yade on the recommendation of my tutor, and compared with PFC, Yade can perform parallel calculation and the calculation speed should be faster, so I asked if there was something wrong with my code.
Thank you1

Revision history for this message
Luc Scholtès (luc) said :
#11

Hello Zhicheng Gao,

if you want to simulate uniaxial compression of rock samples, I'd suggest you have a look at previous related threads/questions.

Please see for example:

https://answers.launchpad.net/yade/+question/661250

https://answers.launchpad.net/yade/+question/694253

https://answers.launchpad.net/yade/+question/635871

Can you help with this problem?

Provide an answer of your own, or ask Zhicheng Gao for more information if necessary.

To post a message you must log in.