yade-users team mailing list archive

[Zhicheng Gao <question699362@xxxxxxxxxxxxxxxxxxxxx>]

Question #699362 on Yade changed:
https://answers.launchpad.net/yade/+question/699362

    Status: Answered => Open

Zhicheng Gao is still having a problem:
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')

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.