New question #669006 on Yade:
https://answers.launchpad.net/yade/+question/669006

Hi,
Recently I wrote a code for traxial simulation. When I try to calculate the 
confining stress, I wrote a subroutine, when I use the Pyrunner(command()), it 
can't work. My whole code is flowing, The problem can be seen in StepNum 1 , 
the Wallstress still be 0 as it was set to be 0 at the beginning :



#################################################
#################################################
#################################################
#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)
#gravel.strength = 200.0e0 # Pa crushable 1000MPa
#gravel.IsSplitable = True
#The rock real young 5.5E10 (5.5e8 also acceptable)

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'

#O.bodies.append(utils.facet(vertices=((-0.6,-0.6,0),(0.6,-0.6,0),(0.6,0.6,0)),dynamic=None,fixed=True,wire=True,color=(0.15,0.15,0.15),highlight=False,noBound=False,material=steel,mask=1,chain=-1))
#O.bodies.append(utils.facet(vertices=((-0.6,-0.6,0),(-0.6,0.6,0),(0.6,0.6,0)),dynamic=None,fixed=True,wire=True,color=(0.15,0.15,0.15),highlight=False,noBound=False,material=steel,mask=1,chain=-1))

### around wall// 30cm diameter// 60cm height
#1#(0.15*0.866,0.15*0.5,z),(0.15*0.5,0.15*0.866,z)//id range (1,10)
#2#(0.15*0.866,0.15*0.5,z),(0.15*0.866,-0.15*0.5,z)//id range (11,20)
#3#(0.15*0.866,-0.15*0.5,z),(0.15*0.5,-0.15*0.866,z)//id range (21,30)
##range is left close, right open

###########make four walls to test the results  ==> each wall height is 0.1

###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
xratio=1
yratio=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'),
]

##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
        ######
        #Step1=> add the loadingplate
        #Step2=> apply the initial axial force and confing force
        #Step3=> apply the loading force and confining stress
        if StepNum == 1:
                PyRunner(command='WallStressGet()',iterPeriod=1)
                #checker.command='WallStressGet()' #get the wall stress
                print "1-Last ballast pos= %.5f,unbalanced forces = %.5f,Wall 
Stress= 
%.5f"%(O.bodies[NumEndBall].state.pos[2],utils.unbalancedForce(),WallStress)
                if O.iter < 30000: return
                if utils.unbalancedForce() > 0.1: 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 == 2:      
                LoadPos=O.bodies[NumLoad].state.pos[2]
                PyRunner(command='WallStressGet()',iterPeriod=1)
                print "2-Loadplate force= %.5f, Wall Stress= 
%.5f"%(plateF,WallStress)
                
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='xyXYZ'
                        StepNum=StepNum+1
                #O.pause()
        elif StepNum == 3:                      
                PyRunner(command='WallStressGet()',iterPeriod=1)
                ConfDevi=(abs(WallStress-ConfStress))/ConfStress ### Unit/kPa # 
parameter!!!
                if ConfDevi > 0.05:
                        
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                
                        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='xyXYZ'
                        MoveVel=WallStress-ConfStress #>0
                        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=(0.000005*MoveVel*0.107*xratio,0.000005*MoveVel*0.107*yratio,0)
                        for i in range(1,57):
                                O.bodies[i].state.blockedDOFs='xyzXYZ'          
        
                        print "6-Caliniconfining stress= %.2f, ConfDevi= %.5f, 
Velocity= %.8f"%(WallStress, ConfDevi, O.bodies[1].state.vel[1])
                else:
                        for i in range(1,57):
                                O.bodies[i].state.blockedDOFs='xyzXYZ'
                        StepNum=StepNum+1
                        IniTime=O.time
                        print "Initial confining stress= %.2f, Step= 
%.2f"%(WallStress, StepNum)
        elif StepNum == 4:
                IniLoadPos=LoadPos
                StepNum=StepNum+1
        elif StepNum == 5:
                print "7-Loadplate pos= %.5f, Loadplate force= 
%.5f"%(O.bodies[NumLoad].state.pos[2], plateF)
                
O.engines=O.engines+[PyRunner(command='Confining()',iterPeriod=1)]+[PyRunner(command='AxialLoading()',iterPeriod=1)]+[PyRunner(command='addPlotData()',iterPeriod=1)]

##ConfiningWall=>wall(1-5)//wall(2-6)//wall(3-7)//wall(4-8)
def WallStressGet():
        global Wall1Stressx,Wall2Stressx,Wall3Stressx,Wall4Stressx
        global Wall1Stressy,Wall2Stressy,Wall3Stressy,Wall4Stressy
        global Wall1S,Wall2S,Wall3S,Wall4S
        global WallStress,ConfStress,ConfDevi,MoveVel,ConfDevi #stress control
        global A1,A2,A3,A4
        global LoadPos,NumLoad,NumEndBall
        ##get the forces
        #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)

#float not ^ just **
        ############## keep confining pressure

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
        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:
                MoveVel=WallStress-ConfStress
                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=(0.000005*MoveVel*0.107*xratio,0.000005*MoveVel*0.107*yratio,0)
                for i in range(1,57):
                        O.bodies[i].state.blockedDOFs='xyzXYZ'
                print "7-Calbrate-loading confining stress= %.5f"%(WallStress)
        else:
                for i in range(1,57):
                        O.bodies[i].state.blockedDOFs='xyzXYZ'
                print "7-loading confining stress= %.5f"%(WallStress)
        ############## keep confining pressure
##AxialLoading
def AxialLoading():
        global forceA,plateF,LoadPos,AxiDevi,IniTime,NumLoad,MoveAxial,AreaPlate
        
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=(plateF-forceA)/forceA
        forceA=250+50*sin((20*pi)*(O.time-IniTime)) #10Hz Omega=2*pi*f// 
120+80*sin(Omega*t)//40-200kPa
        if AxiDevi > 0.03:
                print "clabrate-Loadplate pos= %.5f, Loadplate force= 
%.5f"%(O.bodies[NumLoad].state.pos[2], plateF)
                MoveAxial=plateF-forceA #Unit (kPa)//0.0615 is Area of 
loadingplate
                O.bodies[NumLoad].state.vel=(0,0,0.000005*MoveAxial)
                O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
        else:
                print "finial-Loading-Loadplate pos= %.5f, Loadplate force= 
%.5f"%(O.bodies[NumLoad].state.pos[2], plateF)              
                O.bodies[NumLoad].state.blockedDOFs='xyzXYZ'
##Record
def addPlotData():
        global LoadPos,IniLoadPos,NumLoad,forceA
        global theta,thega,WallStress,Vol,AreaPlate
        theta=forceA
        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,T=O.time,Conf=WallStress,TimeLast=O.time,Volume=Vol)

plot.plots={'Thega':('Theta',),'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)
##############################################
#################################################
Thanks!


-- 
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     : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to