New question #669404 on Yade:
https://answers.launchpad.net/yade/+question/669404
Hello,
Following the last question, I found that a bug for using PyRunner[].
I found that at the begining of
O.engines=[...PyRunner[1...],PyRunner[2...,label=2.],], the
PyRunner[2...,label=2.] was called in the process of PyRunner[1...] by using
'label[2].dead=False', the the PyRunner[2...] works continuously,
but if the PyRunner[2...] was called in the process of PyRunner[1...] by using
the O.engines=O.engines+ PyRunner[2...], the PyRunner[2...] works
discontinuously.
is it a bug for using the PyRunner[]?
The following is two use of PyRunner.
###########################################
####1st using of '----O.engines=O.engines+PyRunner[]'---###
#This simulation for triaxial experiment of ballast which size betweeen
30cm~45cm
#Friction angle for 48 degree
from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,timing
from yade import *
import numpy
from pprint import pprint
import random
from random import uniform
#from random import randint
import math
from math import *
global gravel,steel
gravel = FrictMat()
gravel.density = 2600 #kg/m^3
gravel.young = 2e9
gravel.poisson = 0.21 # real 0.21
gravel.frictionAngle = 0.83 #rad radians(48) // change for rad math.radians(31)
steel = FrictMat()
steel.density = 7850 #kg/m^3
steel.young = 10*gravel.young #inital steel was 10*gravel.young
steel.poisson = 0.3
steel.frictionAngle = 0.55 #rad radians(31)
##
bottom_wall=utils.wall(0.00,axis=2,sense=1,material=steel)
O.bodies.append(bottom_wall)
bottom_wall.state.blockedDOFs='xyzXYZ'
###Number for 7 walls
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-1)),(0.3,-0.15,0.1*(i-1)),(0.3,-0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-1)),(-0.3,-0.15,0.1*i),(0.3,-0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)),(0.3,0.15,0.1*(i-1)),(0.3,0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)),(-0.3,0.15,0.1*i),(0.3,0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-1)),(-0.15,0.3,0.1*(i-1)),(-0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-1)),(-0.15,-0.3,0.1*i),(-0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)),(0.15,0.3,0.1*(i-1)),(0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)),(0.15,-0.3,0.1*i),(0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
global Wall1S,Wall2S,Wall3S,Wall4S
Wall1Stressx=0
Wall2Stressx=0
Wall3Stressx=0
Wall4Stressx=0
Wall1Stressy=0
Wall2Stressy=0
Wall3Stressy=0
Wall4Stressy=0
Wall1S=0
Wall2S=0
Wall3S=0
Wall4S=0
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial #stress control
# Area of the confining Wall
global A1,A2,A3,A4
global LoadPos,IniLoadPos,plateF,IniTime,forceA
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial,AreaPlate
#unit:m^2
IniTime=0
plateF=0 #Unit:kPa
LoadPos=0.6
IniLoadPos=LoadPos # (link to Area of Walls)
forceA=200 # Unit:kPa,P=N/A;N=P*0.0615*1000;A=0.0615
A1=LoadPos*0.3
A2=LoadPos*0.3
A3=LoadPos*0.3
A4=LoadPos*0.3
WallStress=0 # Unit:kPa
ConfStress=100 # Unit:kPa
ConfDevi=0
AxiDevi=0
MoveVel=0
MoveAxial=0
AreaPlate=0.09
# Id of different substances
global
NumLoad,NumEndBall,StepNum,NumEnd,xratio,yratio,zratio,NumContact,WallContact
NumLoad=1
NumEndBall=1
StepNum=1
NumEnd=1
xratio=1
yratio=1
zratio=1
NumContact=4
WallContact=1
# Position and other parameters record
######################parameters
sp=pack.SpherePack()
sp.makeCloud((-0.15,-0.15,0),(0.15,0.15,0.8),rMean=0.016,rRelFuzz=0.25)
sp.toSimulation(material=gravel)
NumEndBall=O.bodies[-1].id#Mark Sphere
global iternum
iternum=0
#O.dt=1.0e-6 #Check it!
O.dt=8.0e-6 #Check it!
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_L3Geom_FrictPhys_ElPerfPl(),Law2_ScGeom_FrictPhys_CundallStrack(),Law2_L6Geom_FrictPhys_Linear()],
),
NewtonIntegrator(damping=0.6,gravity=(0,0,-9.81),label='newton'),
PyRunner(command='TraiStep()',iterPeriod=1,label='checker'),
PyRunner(command='LoadAxial100kPa()',iterPeriod=1,label='loadkeep100kPa'),
]
##Fullfill the box
def TraiStep():
global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
global Wall1S,Wall2S,Wall3S,Wall4S
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial #stress
control,WallStress,ConfStress,ConfDevi,MoveVel
global A1,A2,A3,A4
global
LoadPos,NumLoad,NumEndBall,IniLoadPos,plateF,IniTime,forceA,StepNum,NumEnd,iternum,AreaPlate
global xratio,yratio,ztario,NumContact,WallContact
######
#Step1=> add the loadingplate
#Step2=> apply the initial axial force and confing force
#Step3=> apply the loading force and confining stress
if StepNum == 1:
loadkeep100kPa.dead=True
#loadkeep200kPa.dead=True
StepNum=StepNum+1
elif StepNum == 2:
print "2-unbalanced forces = %.5f"%(utils.unbalancedForce())
if O.iter < 30000: return
if utils.unbalancedForce() > 0.01: return
iternum=O.iter
m=max([b.state.pos[2]+b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
O.bodies.append(utils.wall(m,axis=2,sense=0,material=steel))
NumLoad=O.bodies[-1].id
NumEnd=O.bodies[-1].id
LoadPos=O.bodies[NumLoad].state.pos[2]
StepNum=StepNum+1
elif StepNum == 3:
LoadPos=O.bodies[NumLoad].state.pos[2]
print "3-Loadplate force= %.5f"%(plateF)
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(AreaPlate*1000)
#P=F/A=F/(0.0615*1000)=F/61.5 Unit:kPa
LoadPos=O.bodies[NumLoad].state.pos[2]
if plateF < 100:
O.bodies[NumLoad].state.vel=(0,0,-0.005) #100mm/s
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
O.bodies[NumLoad].state.vel=(0,0,0)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
StepNum=StepNum+1
#O.pause(),first I got to the 200kPa axial stress, then keep
loading axial stress
elif StepNum == 4:
loadkeep100kPa.dead=False
StepNum=StepNum+1
#O.pause()
elif StepNum == 5:
LoadPos=O.bodies[NumLoad].state.pos[2]
A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
for i in range(1,15):
Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
Wall1S=Wall1Stressy/A1
Wall1Stressx=0
Wall1Stressy=0
for i in range(15,29):
Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
Wall2S=Wall2Stressy/A2
Wall2Stressy=0
for i in range(29,43):
Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
Wall3S=Wall3Stressx/A3
Wall3Stressx=0
#Wall3Stressy=0
for i in range(43,57):
Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
Wall4S=Wall4Stressx/A4
Wall4Stressx=0
##########################
WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ### Unit(kPa)
ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa #
parameter!!!
for i in range(1,57):
NumContact=NumContact+len(O.bodies[NumLoad].intrs())
WallContact=NumContact/4+1
NumContact=4
MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*WallContact*(8.0e-6))
################check the parameter
# print "Ini-conf-stress= %.5f, Vel= %.8f, WallContact= %.1f,
NumContact= %.1f, MoveVel= %.8f, Area= %.5f"%(WallStress,
O.bodies[1].state.vel[1], WallContact, NumContact, MoveVel, A1)
################
#MoveVel=0.000005*(WallStress-ConfStress)
if abs(MoveVel) > 0.0001:
MoveVel=0.000001*(WallStress-ConfStress)
else:
print "MoveVel is OK"
for i in range(1,57):
xratio=(abs(O.bodies[i].state.pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
yratio=(abs(O.bodies[i].state.pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
O.bodies[i].state.vel=(MoveVel*xratio,MoveVel*yratio,0)
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Ini-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
O.bodies[1].state.vel[1])
if ConfDevi > 0.05: return
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Ini-Conf-stress= %.5f"%(WallStress)
StepNum=StepNum+1
elif StepNum == 6:
loadkeep100kPa.dead=True
O.engines=O.engines+[PyRunner(command='Confining()',iterPeriod=1)]
StepNum=StepNum+1
O.pause()
elif StepNum == 7:
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=200 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
zratio=1*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
#alpha=50==>100
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
print "force= %.5f, ForceA= %.5f, Vel= %.8f"%(plateF, forceA,
MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
if AxiDevi > 0.05: return
print "Loadplate force= %.5f, ForceA= %.5f"%(plateF, forceA)
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
IniLoadPos=LoadPos
IniTime=O.time
StepNum=StepNum+1
#O.pause()
elif StepNum == 8:
print "8-force= %.5f"%(plateF)
O.engines=O.engines+[PyRunner(command='AxialLoading()',iterPeriod=1)]+[PyRunner(command='addPlotData()',iterPeriod=1)]
StepNum=StepNum+1
O.pause()
else:
print "Well Done"
#O.pause()
##
def Confining():
global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
global Wall1S,Wall2S,Wall3S,Wall4S
global WallStress,ConfStress,ConfDevi,MoveVel #stress control
global A1,A2,A3,A4
global LoadPos,NumLoad,NumEndBall
global xratio,yratio,NumContact,WallContact #control the velocity of
confining walls
LoadPos=O.bodies[NumLoad].state.pos[2]
A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
for i in range(1,15):
Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
Wall1S=Wall1Stressy/A1
Wall1Stressx=0
Wall1Stressy=0
for i in range(15,29):
Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
Wall2S=Wall2Stressy/A2
Wall2Stressy=0
for i in range(29,43):
Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
Wall3S=Wall3Stressx/A3
Wall3Stressx=0
#Wall3Stressy=0
for i in range(43,57):
Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
Wall4S=Wall4Stressx/A4
Wall4Stressx=0
##########################
WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ### Unit(kPa)
ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa #
parameter!!!
if ConfDevi > 0.05:
for i in range(1,57):
NumContact=NumContact+len(O.bodies[NumLoad].intrs())
WallContact=NumContact/4+1
NumContact=4
MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*WallContact*(8.0e-6))
if abs(MoveVel) > 0.0001:
MoveVel=0.000001*(WallStress-ConfStress)
else:
print "MoveVel is OK"
for i in range(1,57):
xratio=(abs(O.bodies[i].state.pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
yratio=(abs(O.bodies[i].state.pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
O.bodies[i].state.vel=(MoveVel*xratio,MoveVel*yratio,0)
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Keep-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
MoveVel)
else:
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Keep-Conf= %.5f"%(WallStress)
############## keep confining pressure
def LoadAxial100kPa():
global
forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,AreaPlate,zratio
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=100 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
zratio=5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
if AxiDevi > 0.05:
print "keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
%.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))#Unit
(kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
print "Done-keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
%.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
##AxialLoading
def AxialLoading():
global
forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,AreaPlate,zratio
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=200+50*sin((20*pi)*(O.time-IniTime)) #10Hz Omega=2*pi*f//
120+80*sin(Omega*t)//40-200kPa
zratio=0.5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
if AxiDevi > 0.05:
print "final-force= %.5f, ForceA= %.5f, Vel= %.8f"%(plateF,
forceA, MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
print "final-Loadplate force= %.5f, ForceA= %.5f"%(plateF,
forceA)
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
##Record
def addPlotData():
global LoadPos,IniLoadPos,NumLoad,forceA,plateF
global theta,thega,WallStress,Vol,AreaPlate
theta=forceA
theta2=plateF
LoadPos=O.bodies[NumLoad].state.pos[2]
thega=((IniLoadPos-LoadPos)/IniLoadPos)*100
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
Vol=IniLoadPos*AreaPlate-LoadPos*(O.bodies[1].state.pos[1])*(O.bodies[1].state.pos[1])
plot.addData(Thega=thega,Theta=theta,Thega2=thega,Theta2=theta2,T=O.time,Conf=WallStress,TimeLast=O.time,Volume=Vol)
##ConfiningWall=>wall(1-5)//wall(2-6)//wall(3-7)//wall(4-8)
plot.plots={'Thega':('Theta',),'Thega2':('Theta2',),'T':('Conf',),'TimeLast':('Volume',)}
plot.plot()
qt.Controller()
V = qt.View()
V.screenSize = (550,450)
V.sceneRadius = 1
V.eyePosition = (0.7,0.5,0.1)
V.upVector = (0,0,1)
V.lookAt = (0.15,0.15,0.1)
#########################################################
#####----2rd is the using of 'label.dead=True'----#################
#This simulation for triaxial experiment of ballast which size betweeen
30cm~45cm
#Friction angle for 48 degree
from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,timing
from yade import *
import numpy
from pprint import pprint
import random
from random import uniform
#from random import randint
import math
from math import *
##################################
#material:ballast and loadingplate
global gravel
global steel
gravel = FrictMat()
gravel.density = 2600 #kg/m^3
gravel.young = 2e9
gravel.poisson = 0.21 # real 0.21
gravel.frictionAngle = 0.83 #rad radians(48) // change for rad math.radians(31)
steel = FrictMat()
steel.density = 7850 #kg/m^3
steel.young = 10*gravel.young #inital steel was 10*gravel.young
steel.poisson = 0.3
steel.frictionAngle = 0.55 #rad radians(31)
##next
#################################################################
##### make circle dormetory
### bottom wall
bottom_wall=utils.wall(0.00,axis=2,sense=1,material=steel)
O.bodies.append(bottom_wall)
bottom_wall.state.blockedDOFs='xyzXYZ'
###Number for 7 walls
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-1)),(0.3,-0.15,0.1*(i-1)),(0.3,-0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((-0.3,-0.15,0.1*(i-1)),(-0.3,-0.15,0.1*i),(0.3,-0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)),(0.3,0.15,0.1*(i-1)),(0.3,0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((-0.3,0.15,0.1*(i-1)),(-0.3,0.15,0.1*i),(0.3,0.15,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-1)),(-0.15,0.3,0.1*(i-1)),(-0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((-0.15,-0.3,0.1*(i-1)),(-0.15,-0.3,0.1*i),(-0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
for i in range(1,8):
O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)),(0.15,0.3,0.1*(i-1)),(0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
O.bodies.append(utils.facet(vertices=((0.15,-0.3,0.1*(i-1)),(0.15,-0.3,0.1*i),(0.15,0.3,0.1*i)),dynamic=None,fixed=True,wire=True,color=(0.35,0.35,0.35),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
global Wall1S,Wall2S,Wall3S,Wall4S
Wall1Stressx=0
Wall2Stressx=0
Wall3Stressx=0
Wall4Stressx=0
Wall1Stressy=0
Wall2Stressy=0
Wall3Stressy=0
Wall4Stressy=0
Wall1S=0
Wall2S=0
Wall3S=0
Wall4S=0
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial #stress control
# Area of the confining Wall
global A1,A2,A3,A4
global LoadPos,IniLoadPos,plateF,IniTime,forceA
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial,AreaPlate
#unit:m^2
IniTime=0
plateF=0 #Unit:kPa
LoadPos=0.6
IniLoadPos=LoadPos # (link to Area of Walls)
forceA=200 # Unit:kPa,P=N/A;N=P*0.0615*1000;A=0.0615
A1=LoadPos*0.3
A2=LoadPos*0.3
A3=LoadPos*0.3
A4=LoadPos*0.3
WallStress=0 # Unit:kPa
ConfStress=100 # Unit:kPa
ConfDevi=0
AxiDevi=0
MoveVel=0
MoveAxial=0
AreaPlate=0.09
# Id of different substances
global NumLoad,NumEndBall,StepNum,NumEnd
NumLoad=1
NumEndBall=1
StepNum=1
NumEnd=1
global xratio,yratio,zratio,NumContact,WallContact
xratio=1
yratio=1
zratio=1
NumContact=4
WallContact=1
# Position and other parameters record
######################parameters
sp=pack.SpherePack()
sp.makeCloud((-0.15,-0.15,0),(0.15,0.15,0.4),rMean=0.016,rRelFuzz=0.25)
sp.toSimulation(material=gravel)
NumEndBall=O.bodies[-1].id#Mark Sphere
global iternum
iternum=0
#O.dt=1.0e-6 #Check it!
O.dt=8.0e-6 #Check it!
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_L3Geom_FrictPhys_ElPerfPl(),Law2_ScGeom_FrictPhys_CundallStrack(),Law2_L6Geom_FrictPhys_Linear()],
),
NewtonIntegrator(damping=0.6,gravity=(0,0,-9.81),label='newton'),
PyRunner(command='TraiStep()',iterPeriod=1,label='checker'),
PyRunner(command='LoadAxial100kPa()',iterPeriod=1,label='loadkeep100kPa'),
PyRunner(command='AxialLoading()',iterPeriod=1,label='axialload'),
PyRunner(command='addPlotData()',iterPeriod=1,label='plotdata'),
PyRunner(command='Confining()',iterPeriod=1,label='keepconf'),
#PyRunner(command='LoadAxial200kPa()',iterPeriod=1,label='loadkeep200kPa'),
]
##Fullfill the box
def TraiStep():
global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
global Wall1S,Wall2S,Wall3S,Wall4S
global WallStress,ConfStress,ConfDevi,MoveVel,AxiDevi,MoveAxial #stress
control,WallStress,ConfStress,ConfDevi,MoveVel
global A1,A2,A3,A4
global
LoadPos,NumLoad,NumEndBall,IniLoadPos,plateF,IniTime,forceA,StepNum,NumEnd,iternum,AreaPlate
global xratio,yratio,ztario,NumContact,WallContact
######
#Step1=> add the loadingplate
#Step2=> apply the initial axial force and confing force
#Step3=> apply the loading force and confining stress
if StepNum == 1:
loadkeep100kPa.dead=True
axialload.dead=True
plotdata.dead=True
keepconf.dead=True
StepNum=StepNum+1
elif StepNum == 2:
#PyRunner(command='WallStressGet()',iterPeriod=1)
#checker.command='WallStressGet()' #get the wall stress
print "2-unbalanced forces = %.5f"%(utils.unbalancedForce())
if O.iter < 30000: return
if utils.unbalancedForce() > 0.05: return
#O.bodies.append(utils.wall(O.bodies[NumEndBall].state.pos[2]+0.04,axis=2,sense=0,material=steel))
iternum=O.iter
m=max([b.state.pos[2]+b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
O.bodies.append(utils.wall(m,axis=2,sense=0,material=steel))
NumLoad=O.bodies[-1].id
NumEnd=O.bodies[-1].id
LoadPos=O.bodies[NumLoad].state.pos[2]
StepNum=StepNum+1
elif StepNum == 3:
LoadPos=O.bodies[NumLoad].state.pos[2]
#PyRunner(command='WallStressGet()',iterPeriod=1)
print "3-Loadplate force= %.5f"%(plateF)
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(AreaPlate*1000)
#P=F/A=F/(0.0615*1000)=F/61.5 Unit:kPa
LoadPos=O.bodies[NumLoad].state.pos[2]
if plateF < 50:
O.bodies[NumLoad].state.vel=(0,0,-0.005) #100mm/s
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
O.bodies[NumLoad].state.vel=(0,0,0)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
StepNum=StepNum+1
#O.pause(),first I got to the 200kPa axial stress, then keep
loading axial stress
elif StepNum == 4:
loadkeep100kPa.dead=False
StepNum=StepNum+1
#O.pause()
elif StepNum == 5:
LoadPos=O.bodies[NumLoad].state.pos[2]
A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
for i in range(1,15):
Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
Wall1S=Wall1Stressy/A1
Wall1Stressx=0
Wall1Stressy=0
for i in range(15,29):
Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
Wall2S=Wall2Stressy/A2
Wall2Stressy=0
for i in range(29,43):
Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
Wall3S=Wall3Stressx/A3
Wall3Stressx=0
#Wall3Stressy=0
for i in range(43,57):
Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
Wall4S=Wall4Stressx/A4
Wall4Stressx=0
##########################
WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ### Unit(kPa)
ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa #
parameter!!!
for i in range(1,57):
NumContact=NumContact+len(O.bodies[NumLoad].intrs())
WallContact=NumContact/4+1
NumContact=4
MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*WallContact*(8.0e-6))
################check the parameter
if abs(MoveVel) > 0.0001:
MoveVel=0.000001*(WallStress-ConfStress)
else:
print "MoveVel is OK"
for i in range(1,57):
xratio=(abs(O.bodies[i].state.pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
yratio=(abs(O.bodies[i].state.pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
O.bodies[i].state.vel=(MoveVel*xratio,MoveVel*yratio,0)
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Ini-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
O.bodies[1].state.vel[1])
if ConfDevi > 0.05: return
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Ini-Conf-stress= %.5f"%(WallStress)
StepNum=StepNum+1
elif StepNum == 6:
loadkeep100kPa.dead=True
#O.engines=O.engines+[PyRunner(command='Confining()',iterPeriod=1)]
keepconf.dead=False
StepNum=StepNum+1
#O.pause()
elif StepNum == 7:
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=200 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
zratio=1*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
#alpha=50==>100
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
print "force= %.5f, ForceA= %.5f, Vel= %.8f"%(plateF, forceA,
MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
if AxiDevi > 0.05: return
print "Loadplate force= %.5f, ForceA= %.5f"%(plateF, forceA)
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
IniLoadPos=LoadPos
IniTime=O.time
StepNum=StepNum+1
#O.pause()
elif StepNum == 8:
print "8-force= %.5f"%(plateF)
#O.engines=O.engines+[PyRunner(command='AxialLoading()',iterPeriod=1)]+[PyRunner(command='addPlotData()',iterPeriod=1)]
axialload.dead=False
plotdata.dead=False
StepNum=StepNum+1
O.pause()
else:
print "Well Done"
#O.pause()
##
def Confining():
global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
global Wall1S,Wall2S,Wall3S,Wall4S
global WallStress,ConfStress,ConfDevi,MoveVel #stress control
global A1,A2,A3,A4
global LoadPos,NumLoad,NumEndBall
global xratio,yratio,NumContact,WallContact #control the velocity of
confining walls
LoadPos=O.bodies[NumLoad].state.pos[2]
A1=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A2=LoadPos*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
A3=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
A4=LoadPos*(O.bodies[22].state.pos[1]-O.bodies[8].state.pos[1])
for i in range(1,15):
Wall1Stressx=Wall1Stressx+abs(O.forces.f(i)[0])
Wall1Stressy=Wall1Stressy+abs(O.forces.f(i)[1])
Wall1S=Wall1Stressy/A1
Wall1Stressx=0
Wall1Stressy=0
for i in range(15,29):
Wall2Stressy=Wall2Stressy+abs(O.forces.f(i)[1])
Wall2S=Wall2Stressy/A2
Wall2Stressy=0
for i in range(29,43):
Wall3Stressx=Wall3Stressx+abs(O.forces.f(i)[0])
Wall3S=Wall3Stressx/A3
Wall3Stressx=0
#Wall3Stressy=0
for i in range(43,57):
Wall4Stressx=Wall4Stressx+abs(O.forces.f(i)[0])
Wall4S=Wall4Stressx/A4
Wall4Stressx=0
##########################
WallStress=(Wall1S+Wall2S+Wall3S+Wall4S)/4000 ### Unit(kPa)
ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa #
parameter!!!
if ConfDevi > 0.05:
for i in range(1,57):
NumContact=NumContact+len(O.bodies[NumLoad].intrs())
WallContact=NumContact/4+1
NumContact=4
MoveVel=1*A1*(WallStress-ConfStress)/((1.0e9)*WallContact*(8.0e-6))
if abs(MoveVel) > 0.0001:
MoveVel=0.000001*(WallStress-ConfStress)
else:
print "MoveVel is OK"
for i in range(1,57):
#zratio=0.5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
xratio=(abs(O.bodies[i].state.pos[0]))/(O.bodies[i].state.pos[0]+0.00001)
yratio=(abs(O.bodies[i].state.pos[1]))/(O.bodies[i].state.pos[1]+0.00001)
O.bodies[i].state.vel=(MoveVel*xratio,MoveVel*yratio,0)
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Keep-Cal-conf-stress= %.5f, Vel= %.8f"%(WallStress,
MoveVel)
else:
for i in range(1,57):
O.bodies[i].state.blockedDOFs='xyzXYZ'
print "Keep-Conf= %.5f"%(WallStress)
############## keep confining pressure
def LoadAxial100kPa():
global
forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,AreaPlate,zratio
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=100 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
zratio=5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
if AxiDevi > 0.05:
print "keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
%.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))#Unit
(kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
print "Done-keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
%.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
def LoadAxial200kPa():
global
forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,AreaPlate,zratio
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=200 #10Hz Omega=2*pi*f// 120+80*sin(Omega*t)//40-200kPa
zratio=5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
if AxiDevi > 0.05:
print "keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
%.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))#Unit
(kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
print "Done-keep-force= %.5f, ForceA= %.5f, Vel= %.8f, CNum=
%.1f"%(plateF, forceA, MoveAxial, len(O.bodies[NumLoad].intrs()))
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
##AxialLoading
def AxialLoading():
global
forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,AreaPlate,zratio
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
plateF=(O.forces.f(NumLoad)[2])/(1000*AreaPlate)
LoadPos=O.bodies[NumLoad].state.pos[2]
AxiDevi=abs((plateF-forceA))/forceA
forceA=200+50*sin((20*pi)*(O.time-IniTime)) #10Hz Omega=2*pi*f//
120+80*sin(Omega*t)//40-200kPa
zratio=0.5*AreaPlate/((2.0e9)*(len(O.bodies[NumLoad].intrs())+1)*(8.0e-6))
MoveAxial=1*zratio*(plateF-forceA)
if abs(MoveAxial) > 0.0001:
MoveAxial=0.000001*(plateF-forceA)
else:
print "MoveAxial is OK"
if AxiDevi > 0.05:
print "final-force= %.5f, ForceA= %.5f, Vel= %.8f"%(plateF,
forceA, MoveAxial)#Unit (kPa)//0.0615 is Area of loadingplate
O.bodies[NumLoad].state.vel=(0,0,MoveAxial)
O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
else:
print "final-Loadplate force= %.5f, ForceA= %.5f"%(plateF,
forceA)
O.bodies[NumLoad].state.blockedDOFs='xyXYZ'
##Record
def addPlotData():
global LoadPos,IniLoadPos,NumLoad,forceA,plateF
global theta,thega,WallStress,Vol,AreaPlate
theta=forceA
theta2=plateF
LoadPos=O.bodies[NumLoad].state.pos[2]
thega=((IniLoadPos-LoadPos)/IniLoadPos)*100
AreaPlate=(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])*(O.bodies[50].state.pos[0]-O.bodies[36].state.pos[0])
Vol=IniLoadPos*AreaPlate-LoadPos*(O.bodies[1].state.pos[1])*(O.bodies[1].state.pos[1])
plot.addData(Thega=thega,Theta=theta,Thega2=thega,Theta2=theta2,T=O.time,Conf=WallStress,TimeLast=O.time,Volume=Vol)
##ConfiningWall=>wall(1-5)//wall(2-6)//wall(3-7)//wall(4-8)
plot.plots={'Thega':('Theta',),'Thega2':('Theta2',),'T':('Conf',),'TimeLast':('Volume',)}
plot.plot()
qt.Controller()
V = qt.View()
V.screenSize = (550,450)
V.sceneRadius = 1
V.eyePosition = (0.7,0.5,0.1)
V.upVector = (0,0,1)
V.lookAt = (0.15,0.15,0.1)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
