yade-users team mailing list archive

Thread
Date
[Question #669006]: Pyrunner(command) problem for Triaxial simulation

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: De zhang <question669006@xxxxxxxxxxxxxxxxxxxxx>
Date: Mon, 14 May 2018 12:52:42 -0000
Reply-to: question669006@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
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.