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