# -*- coding: utf-8 -*- from yade import pack from yade import utils import math import os #from utils import * num_spheres= 10000 mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.1,0.1,0.1) # Top a bottom point of the cube thick = 0.01 compFricDegree = table.compFricDegree finalFricDegree=35 #targetPorosity=0.382 rate=0.01 damp=0.06 stabilityThreshold=1e-3 # leave it too low can have a super stable state, but will eat onion for waiting so long key='_sample' # just used for adding name to the output, kind of boring ## create material #0, which will be used as default O.materials.append(FrictMat(young=356e6,poisson=0.42,frictionAngle=radians(compFricDegree),density=3000,label='spheres')) # hat O.materials.append(FrictMat(young=356e6,poisson=0.5,frictionAngle=0,density=0,label='walls')) # boundary ## create walls around the packing walls=utils.aabbWalls([mn,mx],thickness=thick,oversizeFactor=3,material='walls') wallIds=O.bodies.append(walls) sp=pack.SpherePack() #psdSizes=[0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.0095] # (sizes or radii of the grains vary from 2mm to 9.5mm) #psdCumm=[0.01,0.09,0.25,0.50,0.69,0.90,0.95,1] # the correspondent amount (percentage) of each diameter #psdSizes,psdCumm=[0.001978,0.00218,0.00249,0.002744,0.002894,0.002999,0.003162,0.003305,0.003455,0.00358,0.003709,0.003911,0.004016,0.004273,0.004427,0.004669,0.004968,0.005193,0.005476,0.005826,0.006036,0.00631,0.006595,0.006833,0.007079,0.007335,0.007532,0.007943,0.008526,0.0092],[0.001,0.002,0.005,0.026385,0.047493,0.068602,0.092348,0.12372,0.150396,0.174142,0.200528,0.242744,0.261214,0.311346,0.356201,0.401055,0.469657,0.522427,0.583113,0.654354,0.71504,0.770449,0.823219,0.878628,0.926121,0.968338,0.98153,0.992084,0.997361,1] #--------------------------------------------- #sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,psdSizes,psdCumm,False,seed=1) if num_spheres==1000: sp.makeCloud(mn,mx,0.005,0.65,num_spheres,False,0.8,seed=1) #initial 0.001 for 1000 particle use 0.005 elif num_spheres==10000: sp.makeCloud(mn,mx,0.001,0.65,num_spheres,False,0.8,seed=1) #sp.makeCloud(mn,mx,-1,0,-1,False, 0.95,psdSizes,psdCumm,False,seed=1) sp.toSimulation(material='spheres') volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2]) mean_rad = pow(0.09*volume/num_spheres,0.3333) #clumps=False #if clumps: # c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)]) # sp.makeClumpCloud((-0.24,0.24,-0.24),(0.24,0.24,0.24),[c1],periodic=False) # O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp]) # standalone,clumps=sp.getClumps() # for clump in clumps: # O.bodies.clump(clump) # for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color # #sp.toSimulation() #else: # O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp]) O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away O.usesTimeStepper=True triax=TriaxialStressController( maxMultiplier=1.001, finalMaxMultiplier=1.0001, thickness = thick, radiusControlInterval=10 # stressMask = 7, # max_vel=0.01, # strainDamping=0.99, #0.99 # stressDamping=0.25 ) newton=NewtonIntegrator(damping=damp) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack(label='law')] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()), triax, TriaxialStateRecorder(iterPeriod=50,file='WallStresses'+key,addIterNum=False,initRun=False), newton ] #timestep update interval is 100 by default #O.save('initial'+key+'.xml') #Display spheres with 2 colors for seeing rotations better #Gl1_Sphere.stripes=0 #if nRead==0: yade.qt.Controller(), yade.qt.View() print 'Number of elements: ', len(O.bodies) print 'Box Volume: ', triax.boxVolume print 'Box Volume calculated: ', volume #Display spheres with 2 colors for seeing rotations better Gl1_Sphere.stripes=1 #yade.qt.Controller(), yade.qt.View() ######################################### # Phase 1: ISOTROPIC GENERATOR OF 20kPa # ######################################### while 1: triax.internalCompaction=True # Growing the particles is what we use in this step setContactFriction(radians(compFricDegree)) triax.stressmask=7 triax.goal1=10000 triax.goal2=10000 triax.goal3=10000 O.run(10,True) #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium unb=unbalancedForce() meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3 print 'unbalanced force:',unb,' mean stress: ',meanS # print 'void ratio=',triax.porosity/(1-triax.porosity), 'porosity=', triax.porosity print 'mean stress engine', triax.meanStress, 'kineticE=', utils.kineticEnergy() print '-----State_01: Isotropic compression 10kPa-----' if unb