the bugs in my scripts of layer sample pareparing code
Hi All:
Right now i am doing some simulation based on the following: (i am doing the same job as Nima Goudarzi mentioned)
As a part of a trixal test with my own implemented model in YADE, I need to compact my sample in the vertical direction to reach to a target void ratio. Dimensions of the sample are 3.5*3.5mm*7.1mm and it will be confined in six frictionless walls. The method of compaction is UCM and the sample will be compacted in 5 layers. As I am going to produce 100000 particles (based on a PSD), I may assume that each layer has 2000 spheres in it. Here is the sequence of compacting:
1- First layer is deposited and then compacted to a target void ratio of 0.73.
2- After depositing the second layer on top of the first layer, both layers will be compacted to reach a target void ratio of 0.71.
3- After depositing the third layer on top of the second layer, all three layers will be compacted to reach a target void ratio of 0.69.
4- After depositing the fourth layer on top of the third layer, all four layers will be compacted to reach a target void ratio of 0.67.
5- After depositing the fifth layer on top of the fourth layer, all five layers will be compacted to reach a target void ratio of 0.65.
This 0.65 is the target void ratio for the whole assembly.
the following is the code:
# -*- coding: utf-8 -*-
#######
### DEFINING VARIABLES AND MATERIALS ###
#######
# The following 5 lines will be used later for batch execution
nRead=readParam
num_spheres=2000,# number of spheres
compFricDegree = 30, # contact friction during the confining phase
key='_
unknownOk=True
)
#######
from yade.params import table
a = [1,2,3,4,5]
b = [1]
c = [5]
psdSizes,
for t in a:
num_spheres=
key=table.key
compFricDegree = table.compFricD
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.02 # loading rate (strain rate)
damp=0.7 # damping coefficient
stabilityThres
young=100000000 # contact stiffness
mn,mx=
targetPorosity = (0.73-0.
## create materials for spheres and plates
O.materials.
O.materials.
## create walls around the packing
walls=
wallIds=
while t > b:
O.bodies.erase(4) #delect the rigid plate in the bottom
## use a SpherePack object to generate a random loose particles packing
sp=pack.
sp.makeCloud(
O.bodies.
#or alternatively (higher level function doing exactly the same):
#sp.toSimulati
######
### DEFINING ENGINES ###
######
triax=
## TriaxialStressC
maxMultiplier=
finalMaxMultip
thickness = 0,
## switch stress/strain control using a bitmask.
## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
## 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,
internalCompac
)
newton=
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
GlobalStiffnes
triax,
TriaxialStateR
newton
]
#Display spheres with 2 colors for seeing rotations better
#Gl1_Sphere.
#if nRead==0: yade.qt.
######
### APPLYING Z direction PRESSURE for porosity ##
######
#the value of (isotropic) confining stress defines the target stress to be applied in Z directions
triax.goal2=-50000
while 1:
O.run(1000, True)
##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
unb=unbalance
print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
if unb<stabilityTh
break
#O.save(
print "### Compression state saved ###"
######
### REACHING A SPECIFIED POROSITY PRECISELY ###
######
### We will reach a prescribed value of porosity with the REFD algorithm
import sys #this is only for the flush() below
while triax.porosity>
## we decrease friction value and apply it to all the bodies and contacts
compFricDegree = 0.85*compFricDegree
setContactFric
print "\r Friction: ",compFricDegree," 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)
while t < c:
O.bodies.erase(5) #delect the wall in the positive direction of Z
O.save(
print (targetPorosity)
print "### Compacted state saved ###"
O.bodies.
#######
### APPLYING CONFINING PRESSURE for triaxial test ####
#######
#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=
while 1:
O.run(1000, True)
##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
unb=unbalance
print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
if unb<stabilityTh
break
O.save(
print "### Isotropic state saved ###"
#######
### DEVIATORIC LOADING ###
#######
##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalC
## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFrict
##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5 #mask = x*1 + y*2 + z*4, when x and z, 101
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 50kPa as for the initial confinement.
triax.goal1=-50000
triax.goal3=-50000
##we can change damping here.
newton.damping=0.7
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()
#######
### record and plot data ###
#######
from yade import plot
### a function saving variables
def history():
plot.
s11=
s22=
s33=
sinPhi=
#
i=O.iter)
if 1:
## include a periodic engine calling that function in the simulation loop
O.engines=
##O.engines.
else:
## With the line above, we are recording some variables twice. We could in fact replace the previous
## TriaxialRecorder
## by our periodic engine. Uncomment the following line:
O.engines[
O.run(100,True)
### declare what is to plot. "None" is for separating y and y2 axis
#plot.plots=
#plot.plots=
#plot.plots=
plot.plots=
### the traditional triaxial curves would be more like this:
##plot.
## display on the screen (doesn't work on VMware image it seems)
plot.plot()
##### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) #####
## In that case we can still save the data to a text file at the the end of the simulation, with:
plot.saveDataTx
##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.
plot.saveGnuplo
The major problem is when the code finished the first loop, it will stop at the part of reaching target porosity. The following simulation cannot go through. However, i cannot find the reason. if anyone can point out the bug part, please reply me! Thank you so much and sorry for taking your time!
Regards,
Ting
Question information
- Language:
- English Edit question
- Status:
- Expired
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Last query:
- Last reply: