# coding: utf-8
# loose specimen
from yade import pack
import os
#from utils import *
nRead=utils.readParamsFromTable(
num_spheres=10000,# number of spheres
compFricDegree = 32, # contact friction during the confining phase the final porosity
unknownOk=True,
isoForce=20000,
conStress=100000
)
from yade.params import table
num_spheres=table.num_spheres
mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.1,0.1,0.1)
thick = 0.1
compFricDegree = table.compFricDegree
finalFricDegree=35
#targetPorosity=0.382
rate=0.05
damp=0.06
stabilityThreshold=0.001 # leave it too low can have a super stable state, but will eat onion for waiting so long
#stopStep=160000 # pas de calcul on choisit pour faire calculer travail du second-ordre
key='_probing_envelope_' # 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'))
O.materials.append(FrictMat(young=356e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))
## create walls around the packing
walls=utils.aabbWalls([mn,mx],thickness=thick,oversizeFactor=5,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
#sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,psdSizes,psdCumm,False,seed=1)
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,
)
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), # 20 is too small, will make big file size
newton
]
#O.save('initial'+key+'.xml')
#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=1
#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
# Applying confining force
#########################################
# 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=table.isoForce
triax.goal2=table.isoForce
triax.goal3=table.isoForce
O.run(1000,True)
#the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
unb=unbalancedForce()
#average stress
#note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
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 20kPa--(^_^)---'
if unb<stabilityThreshold and abs(meanS-table.isoForce)/table.isoForce<0.001:
break
O.save('confinedState'+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)
###### porosity reaching #####
#import sys #this is only for the flush() below
#while triax.porosity>targetPorosity:
# # we decrease friction value and apply it to all the bodies and contacts
# compFricDegree = 0.95*compFricDegree
# setContactFriction(radians(compFricDegree))
# print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
# sys.stdout.flush()
# # while we run steps, triax will tend to grow particles as the packing
# # keeps shrinking as a consequence of decreasing friction. Consequently
# # porosity will decrease
# O.run(500,1)
#O.save('compactedState'+key+'.yade.gz')
#print "### Compacted state saved ###"
################################################################################
# 2nd phase: Confining the specimen to desired value of stress of confinement,
# reduce the cinetic energy
# will implement in this state the code for generating the granulometric curve
################################################################################
while 1:
triax.internalCompaction=False
triax.stressMask=7
setContactFriction=radians(compFricDegree)
triax.goal1=table.conStress
triax.goal2=table.conStress
triax.goal3=table.conStress
O.run(1000, True)
#the global unbalanced force on dynamic bodies, thus excluding boundaries,
#which are not at equilibrium
unb=unbalancedForce()
#average stress
#note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
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 '-----Current State_02: Compress 100kPa--v(-___-)v--'
if unb<stabilityThreshold and abs(meanS-table.conStress)/table.conStress<0.001:
break
O.save('compressState'+key+'.xml') # important, put to meromy this file to start the deviatoric loading later
print 'step that starts the deviatoric loading ', O.iter
e22Check=triax.strain[1]
# We will save a check point, to know that from this point the deviatoric loading process begins, so we know where to start the ploting shit
from yade import plot
from pprint import pprint
plot.reset()
plot.addData(CheckpointStep=O.iter,Porosity=triax.porosity,e11=triax.strain[0],e22=triax.strain[1],e33=triax.strain[2],elasticE=law.elasticEnergy(),kineticE=utils.kineticEnergy())
plot.saveDataTxt('checkpoint.txt',vars=('CheckpointStep','Porosity','e11','e22','e33','elasticE','kineticE'))
###############################################
# 3rd phase: DEVIATORIC LOADING (BIG SHOW!!!) #
###############################################
#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction=radians(35)
#TriaxialStateRecorder(iterPeriod=20,file='WallStresses_load_deviatoric'+key,addIterNum=True,initRun=True)
while 1:
triax.stressMask=5
triax.goal2=-rate
triax.goal1=triax.goal3=table.conStress
O.run(100, True)
#the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
unb=unbalancedForce()
#note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
axialS=triax.stress(triax.wall_top_id)[1]
print 'step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-table.conStress
print 'axial deformation (%)', (triax.strain[1]-e22Check)*100
if abs((triax.strain[1]-e22Check)-0.01)<=0.001:
O.save('firstpoint.xml')
if abs((triax.strain[1]-e22Check)-0.02)<=0.001:
O.save('secondpoint.xml')
if abs((triax.strain[1]-e22Check)-0.03)<=0.001:
O.save('thirdpoint.xml')
if abs((triax.strain[1]-e22Check)-0.05)<=0.001:
O.save('fourthpoint.xml')
if triax.strain[1]-e22Check>=0.3 : #triax.sigma2
break
O.save('final.xml')
Comments