# -*- 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<stabilityThreshold and abs(meanS-table.isoForce)/table.isoForce<0.001:
break
O.save('confinedState20'+key+'.xml') # remember this impotant part, will load is later
print "## Isotropic state saved (20 kPa) ##"
print 'current porosity=',triax.porosity
print 'current void ratio=',triax.porosity/(1-triax.porosity)
Comments