yade-users team mailing list archive

[wangxiaoliang <question216109@xxxxxxxxxxxxxxxxxxxxx>]

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

wangxiaoliang gave more information on the question:
Well, our mission is to understand the characteristic of granular of
granular media after erosion by water, background of which is

piping or dam-break.  In my opinion, erosion is a lose of small size
particles, so I just want to delete some small particles in a triaxial

compression tests. scripts is below. script will run to erosion
successfully, but after that error was met "Segmentation fault (core
dumped)"

#####################################333
# encoding: utf-8
# a script to simulate the degradation duo to erosion
from yade import pack

# The following 5 lines will be used later for batch execution
nRead=utils.readParamsFromTable(
        num_spheres=500,# number of spheres
        compFricDegree = 30, # contact friction during the confining phase
        unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below)
key='_phi='+str(compFricDegree) # put you simulation's name here
young=100e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,2,1) # corners of the initial packing
thick = 0.01 # thickness of the plates

sp=pack.SpherePack()

## box between mn and mx, avg radius ± ½(20%), 2k spheres
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=-1,rRelFuzz=.2,num=num_spheres,periodic=0,porosity=0.95)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=young,poisson=.4,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=.4,frictionAngle=0,density=0,label='frictionless'))

volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
mean_rad = pow(0.09*volume/num_spheres,0.3333)
mean_rad = 5*mean_rad

#define psd and spheres cloud
psdSizes = [1.*mean_rad,1.81*mean_rad,2.58*mean_rad,3.35*mean_rad,4.13*mean_rad,4.90*mean_rad,5.68*mean_rad,6.45*mean_rad,7.22*mean_rad,8.0*mean_rad]
psdCumm =[0.,.19,.24,.27,.31,.34,.38,.42,.49,1.]
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=num_spheres)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

## create material #0, which will be used as default
## it seems that wall should be created after the spheres
walls=utils.aabbWalls([mn,mx],thickness=1e-10,material='frictionless')
wallIds=O.bodies.append(walls)
#print out psd generated
dmj = utils.psd(bins=9)[0]
hmj = utils.psd(bins=9)[1]
fileHandle = open ( 'psdGenerated'+str(compFricDegree), 'w' )
fileHandle.write('r'+' '+'h'+'\n')
for i in [0,1,2,3,4,5,6,7,8,9]:
        r = dmj[i]
        h = hmj[i]
        fileHandle.write(str(r)+'  '+str(h)+'\n')
fileHandle.close()
# print the psd required in to file
fileHandle = open ( 'psdRequired'+str(compFricDegree), 'w' )
fileHandle.write('r'+' '+'h'+'\n')
for i in [0,1,2,3,4,5,6,7,8,9]:
        r = psdSizes[i]
        h = psdCumm[i]
        fileHandle.write(str(r)+'  '+str(h)+'\n')
fileHandle.close()

triax=TriaxialCompressionEngine(
        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],
        internalCompaction=1,
        epsilonMax = 0.2,
        sigmaIsoCompaction=50e3,
        sigmaLateralConfinement=50e3,
        max_vel=10,
        strainRate=rate,
        noFiles = 0,
        Key = key,
        label="triax"
)

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
#       TriaxialStateRecorder(iterPeriod = 100, file = 'WallStress'),
        # you can add TriaxialStateRecorder and such here…
        NewtonIntegrator(damping=.4)
]
while 1:
        O.run(1000,True)
        print 'O.iter = ',O.iter
        print 'StrainY = ',triax.strain[1]
        if triax.strain[1] > 1e-4:
                break


triax.setContactProperties(finalFricDegree)


################################################3
### start to erosion
### Now only a test
################################################
print "start erosion\n"
radius_erosion = dmj[1]
bodiesToBeDeleted = []
#for b in O.bodies:
#       if b.id < 500:
#               if b.shape.radius < radius_erosion:
#                       bodiesToBeDeleted.append(b)
#for b in bodiesToBeDeleted:
#       O.bodies.erase(b.id)
O.bodies.erase(100)
print "Erosion finished...."

###psd after erosion
# print out psd after erosion
dmj = utils.psd(bins=9)[0]
hmj = utils.psd(bins=9)[1]
fileHandle = open ( 'psdAfterErosion'+str(compFricDegree), 'w' )
fileHandle.write('r'+' '+'h'+'\n')
for i in [0,1,2,3,4,5,6,7,8,9]:
        r = dmj[i]
        h = hmj[i]
        fileHandle.write(str(r)+'  '+str(h)+'\n')
fileHandle.close()
                                                  
# print the psd required in to file
fileHandle = open ( 'psdAfterErosion'+str(compFricDegree), 'w' )
fileHandle.write('r'+' '+'h'+'\n')


##################################################
#### Define Fabric Tensor and calculattion at Sep.26
##################################################
def CalFabric():
        F = Vector3(0,0,0)
        c_count = 0
        ev = triax.strain[0]+triax.strain[1]+triax.strain[2] # volume strain
        V = 2*(1.0 - ev) # volume now
        for b in O.interactions:
            if (b.geom.penetrationDepth > 0):
                c_count +=1
                n = b.geom.normal
#               nj = b.geom.normal
#               Ri = b.geom.refR1
#               Rj = b.geom.refR2
                F[0] += n[0]*n[0]
                F[1] += n[1]*n[1]
                F[2] += n[2]*n[2]
        F[0] = F[0]/c_count
        F[1] = F[1]/c_count
        F[2] = F[2]/c_count
        return F
from yade import plot
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1000,command='history()',label='recorder')]+O.engines[5:6]
def history():
        F = Vector3(0,0,0)
        F = CalFabric()
        IF1 = F[0]+F[1]+F[2]
        IF2 = F[1] - 0.5*(F[0] + F[1])
        z = utils.avgNumInteractions(skipFree=True)
        plot.addData(e11=triax.strain[0], e22=triax.strain[1], e33=triax.strain[2],
                    s11=-triax.stress(0)[0],
                    s22=-triax.stress(2)[1],
                    s33=-triax.stress(4)[2],
                    IF1 = IF1,
                    IF2 = IF2,
                    Z = z,
                    n = triax.porosity,
                    i=O.iter)
plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}

O.saveTmp()
#plot.plot()
while 1:
        O.run(1000,True)
        if O.iter%1000 ==0:
                print 'iter=',O.iter
                print 'StrainY = ',triax.strain[1]
        if triax.strain[1] > 0.2:
                break

## In that case we can still save the data to a text file at the the end of the simulation, with: 
plot.saveDataTxt('results'+key)
##or even generate a script for gnuplot. Open another terminal and type  "gnuplot plotScriptKEY.gnuplot:
plot.saveGnuplot('plotScript'+key)

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.