# -*- coding: utf-8 -*-
# This code is based on the example initially created by Nejib Hada, Luc Sibille, Bruno Chareyre, this is a fix modified by Nguyen Nho Gia Hien
# Any further questions please put them on launchpad.net/yade and/or email nho-gia-hien.nguyen@insa-lyon.fr (English, Francais, Tieng Viet)
# in this script the method of confining is moving wall
from yade import pack
#from utils import *
nRead=utils.readParamsFromTable(
num_spheres=5000,# number of spheres (first try with lower value for faster testing calculation, minimum value should be 10000 to have good result)
compFricDegree = 4, # contact friction during the confining phase the final porosity is slightly bigger than 0.409, it means we must lower this slightly, maybe 18.5
unknownOk=True,
isoForce=100000, # stand for the stress isotrope state that we define in the compression state
conStress=100000 # stand for the confinement stress, actually these two parameters are the smae since I developped this for the former ThreeDEngine
)
from yade.params import table
num_spheres=table.num_spheres
## corners of the initial packing
mn,mx=Vector3(-0.2,-0.2,-0.2),Vector3(0.2,0.2,0.2)
thick = 0.1 # thickness of the box of the specimen (currently don't see any effect of this value, put BIG if you want to)
compFricDegree = table.compFricDegree
finalFricDegree=35 # the angle of friction we will use for the deviatoric loading, microscopic angle of friction
targetPorosity=0.382
rate=0.05
damp=0.06
stabilityThreshold=0.001
compForce=100000
stickLoad=130000 # niveau de contrainte = 30000/100000 = 0.3
key='_probing_envelope_'
## 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=18e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))
## create walls around the packing
walls=utils.aabbWalls([mn,mx],thickness=thick,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.65,psdSizes,psdCumm,False,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)
## Uncomment this section if we consider creating clump pack instead of separated simple particles
#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, # spheres growing factor (fast growth)
# finalMaxMultiplier=1.00001, # spheres growing factor (slow growth)
# turn the maxMultiplier and finalmaxMultiplier on if we use the internal growing particles for the confining process
thickness = thick,
radiusControlInterval=20, # this impact makeCloud, 10 means 10 values, 20 mean 20 values and so on, this affect the PSD curve dramatically
# switch stress/strain control using a bitmask
# Binary variety: x=1 if stress is controlled on x, else x=0. Same for for y and z, which are 1=true or 0=false.
# Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
# to put it differently, the mask is the integer whose binary representation is xyz, i.e.
# "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
stressMask = 7, # this means we control stress of all the axes x, y and z
max_vel=1, # the maximum speed of the moving wall, can be reduced with intel to lower the cinetic energy
internalCompaction=False, # If true the confining pressure is generated by growing particles
)
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()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
triax,
TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key), # data will be saved every 100 steps, reduce if we like more detail and heavier data to handle.
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) # this should be = number of particles + 6 walls. For this example, this value will be 1006
print 'Box Volume: ', triax.boxVolume
print 'Box Volume calculated: ', volume
# Applying confining force
triax.goal1=table.isoForce
triax.goal2=table.isoForce
triax.goal3=table.isoForce
while 1:
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
if unb<stabilityThreshold and abs(meanS-table.isoForce)/table.isoForce<0.001:
break
O.save('confinedState'+key+'.xml')
print "### Isotropic state saved ###"
print 'current porosity (n)=',triax.porosity
print 'current void ratio (e)=',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+'.xml')
print "### Compacted state saved ###"
print 'current porosity (n) =',triax.porosity
print 'current void ratio (e)=',triax.porosity/(1-triax.porosity)
print 'step that starts the deviatoric loading ', O.iter
## the following line is a temporary method to know from where to handle the data, will be improved later
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])
plot.saveDataTxt('checkpoint.txt',vars=('CheckpointStep','Porosity','e11','e22','e33'))
##############################
### DEVIATORIC LOADING ###
##############################
#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False # always
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(35))
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
#O.saveTmp()
# From this point we start the probing process
#First perform a transverse (horizontal) isotropic compression (or a stress controlled drained triaxial compression) to reach an initial state from where stress probes will be applied
#... need to active stress control in 3 directions
#... choose the value of axial stress where we want to stop the compression
#... it will be stickLoad variety
## from this we use impose a strain rate on the axial axe and keep the confinement stress constant, until the anisotropic state is attained: stickLoad
while 1:
triax.stressMask=7
triax.goal2=stickLoad
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 'unbalanced force:',unb,' sigma2: ',axialS
# if unb<stabilityThreshold and abs(axialS-(stickLoad-1000))/(stickLoad-1000)<0.001: #triax.sigma2
if unb<stabilityThreshold and abs(axialS-stickLoad)/stickLoad<0.001: #triax.sigma2
break
plot.reset()
plot.addData(CheckpointStep2=O.iter,Porosity2=triax.porosity,e11=triax.strain[0],e22=triax.strain[1],e33=triax.strain[2])
plot.saveDataTxt('checkpoint2.txt',vars=('CheckpointStep2','Porosity2','e11','e22','e33'))
# return to the stress control state, maintain the anisotropic situation
#while 1:
# triax.stressmask=7
# triax.goal2=stickLoad
# triax.goal1=triax.goal3=100000
# O.run(100, True)
# unb=unbalancedForce()
# axialS=triax.stress(triax.wall_top_id)[1]
# print 'unbalanced force (anisotropic)', unb, 'sigma2 (anisotropic): ',axialS
# if unb<stabilityThreshold:
# break
#plot.reset()
#plot.addData(CheckpointStep=O.iter,Porosity=triax.porosity,e11=triax.strain[0],e22=triax.strain[1],e33=triax.strain[2])
#plot.saveDataTxt('checkpoint_anisotrope.txt',vars=('CheckpointStep','Porosity','e11','e22','e33'))
O.save('anisotropicState'+key+'.xml') # save the state so that we can reload after for the probing calculation
#Perform stress probes from the anisotropic state
dSnorm = 10000.0 #norm of the stress increment 0.01 kPa
nbProbes = 18 #number of stress directions tested
rampIte = 20 #nb iterations to increase the stress state until the final desired stress value
#LUC: je fixe des nombres d'iterations c'est moins elegant que de chercher explicitement un etat d'equilibre mais ca permet de poursuivre le calcul meme si un etat de contrainte n'est pas correctement atteint pour un stress probe et qu'il est difficile de se stabiliser a cet etat de contrainte (i.e. attendre longtemps...)
stabIte = 5000 #nb iterations to stabilize sample after reaching the final stress value
# an array for saving stress increments and strain responses; arrays are in "numpy" extension
import numpy
probings=numpy.zeros((3,nbProbes))
def increment(dsr=0,dsa=1):
for ite in range(rampIte):# progressivaly increase of stress state
O.run(20, True)
#incrementation of stress state
triax.goal2 = initSa+dsa/rampIte*ite
triax.goal1 = triax.goal3 = initSr+dsr/rampIte*ite
print triax.goal1, triax.goal2, 'stress value'
# fix the stress value for stabilization at the final state
triax.goal2 = initSa+dsa
triax.goal1 = triax.goal3 = initSr+dsr
while 1:
O.run(100, True)
unb=unbalancedForce()
print 'unbalanced force:',unb,' strain: ',triax.strain
O.save('checkpoint.xml') # reload this file to check if we see mistakes
if unb<stabilityThreshold: break
# loop over all the stress directions
for i in range(nbProbes):
# computation of the stress direction of the current stress probe
alphaS = 2*pi/nbProbes*(i-1)
print 'stress probe nb:',i,' stress direction (deg): ',degrees(alphaS)
# computation of the stress increment in the axial direction
dSa = dSnorm*sin(alphaS)
# computation of the stress increment in the radial direction
dSr = dSnorm*cos(alphaS)/sqrt(2.0)
#Load the initial anisotropic state before running a new stress probe
O.load('anisotropicState'+key+'.xml')
triax.stressMask=7
triax.goal2=stickLoad
triax.goal1=triax.goal3=100000
#We redefine the "triax" label, else it would point to inactive engine from previous simulation that is still in memory
triax=O.engines[4]
initSa=triax.goal2 #save of the initial axial stress
initSr=triax.goal1 #save of the initial radial stress
# define the final stress state to be reached
finalSa = initSa+dSa
finalSr = initSr+dSr
#... need to active stress control in 3 directions if not yet done
# triax.stressControl_1=triax.stressControl_2=triax.stressControl_3=True
triax.stressMask=7 # Hien: the previous is 7, that means the stress control is on, but the following code is focus on strain control so I turn this off
# fix a high value of maximum strain rate, the progressive loading will be done by progressively increasing the desired stress state at each iteration
# triax.strainRate1=triax.strainRate2=triax.strainRate3=1000.0
# triax.goal1=triax.goal2=triax.goal3=10 #!!!!!! # initial value is 1000
increment(dSr,dSa)
probings[:,i]=triax.strain
#open a file for writing probing results
a=open('probings'+key,'w')
for i in range(nbProbes): a.write(str(probings[0][i])+' '+str(probings[1][i])+' '+str(probings[2][i])+'\n')
a.close()
#plot
from pylab import *
plot(probings[0,:]*sqrt(2),probings[1,:],'bo--')
ylabel(r'$\epsilon_{22}$')
xlabel(r'$\epsilon_{11} \times \sqrt{2}$')
title('response envelop')
grid()
show()
O.pause()
when modelling in DEM is needed to define which state the simulation is in: quasi-static or dynamic.
and for dynamic system we’ll have 
(J).
: density of the particle i
: velocity of particles
: inertia tensor
Crocodile & Saviour 