yade-users team mailing list archive

[黎犴dada <question699030@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Needs information => Open

黎犴dada gave more information on the question:
Thank you for your answer. This is the code of my three-axis generation.
When creating particles, I want to use the grading curve I already have
to create particles with specific data, but I am confused about how to
modify this code.

# unicode: UTF-8
# For 2D biaxial simulation
# 21/10/2016
# Yade version 2016.06a-24-0557faf~trusty

#########################################
### Defining parameters and variables ###
#########################################

#Material constants 
Density = 3000
FrictionAngle = 35
PoissonRatio = 0.5
Young = 300e6
Damp = 0.5
AvgRadius = 0.0027
N_particles = 25000

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9

#Packing variables
mn = Vector3(0.,0.,0.)
mx = Vector3(1.5,1.5,1.5)
targetPorosity = 0.2
Porosity = 0.

#Confining variables
ConfPress1 = -90000  #pre-compression
ConfPress = -1.0e5

#Loading control
LoadRate = -0.01

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

########################################
#import necessary packages
from yade import pack,plot,os,timing
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

########################################
### Sample Preparing ###################
########################################

#Create materials for spheres and plates
SphereMat = O.materials.append(FrictMat(young = Young, poisson = PoissonRatio, frictionAngle = radians(FrictionAngle), density = Density))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))

#Create walls for packing
wallIds = O.bodies.append(aabbWalls([mn,mx],thickness=0.001,material=WallMat))

#Use SpherePack object to generate a random loose particles packing
#O.periodic = True
#O.cell.setBox(8,3,8)

sp = pack.SpherePack()

#psdSizes,psdCumm=[.003,0.0035,0.004,0.0045,0.005,0.0055,0.006,0.0065,0.007,0.0075,0.008,0.0085,0.009],[0.,0.003,0.025,0.081,0.182,0.325,0.493,0.660,0.8,0.890,0.957,0.984,1.]
#pylab.plot(psdSizes,psdCumm,label='precribed PSD')
sp.makeCloud(Vector3(0,0,0.0),Vector3(1.5,1.5,1.5) ,-1,0.33,N_particles,False, 0.75)
#pylab.plot(*sp.psd(bins=30,mass=True),label='PSD of (free) %d random spheres'%len(sp))

#pylab.legend()

#pylab.show()
sp.toSimulation(material = SphereMat)


O.usesTimeStepper=True
O.trackEnerty=True
#################################
#####Defining triaxil engines####
#################################

###first step: compression#######
triax1=TriaxialStressController(
    #define wall ids
    #wall_bottom_id = wallIds[4],
    #wall_top_id = wallIds[5],
    #wall_left_id = wallIds[1],
    #wall_right_id = wallIds[0],
    #wall_back_id = wallIds[2],
    #wall_front_id = wallIds[3],
    #wall_back_activated = False, #for 2d simulation
    #wall_front_activated = False,
    thickness = 0.001,
    maxMultiplier=1.+1.5e5/Young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+4e3/Young,
    #maxMultiplier = 1.002,
    internalCompaction = True, # If true the confining pressure is generated by growing particles
    #max_vel = 1.5,
    stressMask = 7,
    computeStressStrainInterval = 10,
    
    goal1 = ConfPress1,
    goal2 = ConfPress1,
    goal3 = ConfPress1,
    
)


#O.dt=0.5*PWaveTimeStep()

newton=NewtonIntegrator(damping=Damp)


###engine
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, defaultDt=4*utils.PWaveTimeStep()),
    triax1,
    newton,
    PyRunner(realPeriod=10,command='checkUnbalanced()',label='check'),
    PyRunner(command='addPlotData()',iterPeriod=2000,label='record'),
    TriaxialStateRecorder(iterPeriod=2000,file='WallStresses')
]
# Simulation stop conditions defination
def checkUnbalanced():
    unb=unbalancedForce()
    mStress = (triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1]+triax1.stress(triax1.wall_front_id)[2])/3.
    s1 = triax1.stress(triax1.wall_right_id)[0]
    s2 = triax1.stress(triax1.wall_top_id)[1]
    s3 = triax1.stress(triax1.wall_front_id)[2]
    
    if unb<0.01 and abs(ConfPress1-mStress)/(-ConfPress1)<0.01:
       O.pause()
       O.save('initial.yade.bz2')
       ################################## postprocessing for simulation ######################################################
       f = open("particleinfo.dat",'w')
       f.write('# This is the result data of 2D simulation\n\n')
       f.write('# There are 8 types of varibles in this data as follows:\n\n')
       f.write("varibles='X-cordinate','Y-cordinate','Z-cordinate','Radius','X-displacement','Y-displacement','Z-displacement','Ids'\n\n")
       f.write('%-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s\n'% ('X-cordinate','Y-cordinate','Z-cordinate','Radius','X-displacement','Y-displacement','Z-displacement','Ids','dens'))
       for b in O.bodies:
           if isinstance(b.shape,Sphere):
              pos = b.state.pos
              rad = b.shape.radius
              displ = b.state.displ()
              #vel = b.state.vel #vector3
              #O.forces.f(b) for forces
              #O.forces.t(b) for torques
              f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16d %-16g\n'%(pos[0],pos[1],pos[2],rad,displ[0],displ[1],displ[2],b.id,b.material.density))
       f.write('################################ ends ##########################################')
       f.close()

       
       

# collect history of data
def addPlotData():
    unb = unbalancedForce()
    mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1]+triax1.stress(triax1.wall_front_id)[2])/3.

    
   
    plot.addData(e11=-triax1.strain[0], e22=-triax1.strain[1], e33=-triax1.strain[2],
  		 ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
		 s11=-triax1.stress(triax1.wall_right_id)[0],
		 s22=-triax1.stress(triax1.wall_top_id)[1],
		 s33=-triax1.stress(triax1.wall_front_id)[2],
                 ub=unbalancedForce(),
                 dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
                 disx=triax1.width,
                 disy=triax1.height,
                 disz=triax1.depth,
		 i=O.iter,
                 porosity=utils.porosity(),
                 ke=utils.kineticEnergy(),
                 totalEnergy=O.energy.total()
                 )

    
    print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f,s33:%f, coordination number: %f, porosity: %f' %(unb,mStress,-triax1.stress(triax1.wall_right_id)[0],-triax1.stress(triax1.wall_top_id)[1],-triax1.stress(triax1.wall_front_id)[2],avgNumInteractions(),utils.porosity())) 

    plot.saveDataTxt('loadinglog.txt.bz2')
    #plot.saveGnuplot('plotScript')


###display
#yade.qt.Controller(),yade.qt.View()


### declare what is to plot. "None" is for separating y and y2 axis
plot.plots={'i':('ub',),'i ':('e11','e22','e33',),' i':('s11','s22','s33',),'e22':'dstress'}
#plot.plots={'i':('ub',),'i ':('s11','s22','s33'),' i':('e11','e22','e33')}

### the traditional triaxial curves would be more like this:
##plot.plots={'e22':('s11','s22','s33',None,'ev')}

## display on the screen (doesn't work on VMware image it seems)

#plot.plot()
O.run()#; O.wait()

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