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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp