yade-users team mailing list archive

[Wojciech Sobieski <question255937@xxxxxxxxxxxxxxxxxxxxx>]

Question #255937 on Yade changed:
https://answers.launchpad.net/yade/+question/255937

Wojciech Sobieski posted a new comment:
Hi Bruno!

I sent you wrong script - this is my last try:

#--------------------------------------------------------------------------------------------
from yade import plot
from yade.pack import *

from yade import utils

O=Omega()

#--------------------------------------------------------------------------------------------
#rozmiar obszaru obliczeniowego:
length = 0.25
height = 0.5
width = 0.25

#grubosc scianek (musza to byc bryly a nie plaszczyzny):
thickness = 0.001

#liczba czastek:
number = 5000

#porowatosc docelow:
targetPorosity = 0.41	

#predkosc kompresji (scianki, ktora sie rusza):
v = 5.0

#liczba iteracji:
maxIter = 150000

#plik:
file = 'p06_25'

#sprawdzenie liczby przedzialow:
n_band = sum(1 for line in open(file+'.txt'))

#czytanie danych o rozkladzie zloza:
input = open(file+'.txt','r')
psdSizes = []
psdCumm = []
try:
  for line in input:
    psdSizes.append(float(line[0:16]))    
    psdCumm.append(float(line[18:]))
finally:
  input.close()

#korekta wartosci skrajjnych (musi byc 0 i 1):
psdCumm[0] = 0.0
psdCumm[n_band-1] = 1.0

#informacje o zadaniu:
os.system('clear')
print '' 
print ' length = ', length
print ' height = ', height
print ' width = ', width
print ' thickness = ', thickness
print ' v = ', v
print ' maxIter = ', maxIter
print ' n_bands = ',n_band
print ' psdSizes, psdCumm ' 
print '' 

#dane stali:
stDensity = 7860 #sprawdzone
stYoung = 2e11 #sprawdzone
stPoisson = 0.3 #sprawdzone
stFrAngle = 30 #

#dane szkla:
spDensity = 2600 #sprawdzone
spYoung = 50e9 #sprawdzone
spPoisson = 0.18 #sprawdzone 0.18-0.3
spFrAngle = 5 # 0.449

#--------------------------------------------------------------------------------------------
#definicja materialu scian:
O.materials.append(FrictMat(density=stDensity,young=stYoung,poisson=stPoisson,frictionAngle=radians(stFrAngle),label='wallMat'))

#definicja poszczegolnych bokow kostki tworzacej obszar obliczeniowy:
leftBox = box( center=(-thickness/2.0,(height)/2.0,0), extents=(thickness/2.0,5*(height/2.0+thickness),width/2.0) ,fixed=True,wire=True)
lowBox = box( center=(length/2.0,-thickness/2.0,0), extents=(length/2.0,thickness/2.0,width/2.0) ,fixed=True,wire=True)
rightBox = box( center=(length+thickness/2.0,height/2.0,0), extents=(thickness/2.0,5*(height/2.0+thickness),width/2.0) ,fixed=True,wire=True)
upBox = box( center=(length/2.0,height+thickness/2.0,0), extents=(length/2.0,thickness/2.0,width/2.0) ,fixed=True,wire=True)
behindBox = box( center=(length/2.0,height/2.0,-width/2.0-thickness/2.0), extents=(2.5*length/2.0,height/2.0+thickness,thickness/2.0), fixed=True,wire=True)
inFrontBox = box( center=(length/2.0,height/2.0,width/2.0+thickness/2.0), extents=(2.5*length/2.0,height/2.0+thickness,thickness/2.0), fixed=True,wire=True)

#dodawanie cial do sceny - mozna to uzyc w state.pos[]
O.bodies.append([leftBox,lowBox,rightBox,upBox,behindBox,inFrontBox])

#--------------------------------------------------------------------------------------------
#definicja materialu sfer:
O.materials.append(FrictMat(density=spDensity,young=spYoung,poisson=spPoisson,frictionAngle=radians(spFrAngle),label='spMat'))

#definiowanie pozycji chmury:
mincorner = Vector3(0,0.0,-width/2.0)
maxcorner = Vector3(length,height,width/2.0)

#tworzenie chmury:
sp = yade._packSpheres.SpherePack()
sp.makeCloud(mincorner,maxcorner,psdSizes=psdSizes,num=number,psdCumm=psdCumm,distributeMass=1)

#dodanie chmury do sceny:
O.bodies.append([sphere(s[0],s[1],material='spMat') for s in sp])

#--------------------------------------------------------------------------------------------
#informacja o liczbie sfer i o porowatosci:
print ' particles = ',len(O.bodies)
print ' porosity = ',utils.porosity(length*height*width) 

plot.resetData()

#---------------------------------------------------------------------------------
triax = TriaxialStressController(
  thickness = 0,
  stressMask = 7,
  internalCompaction = True,
  )

#--------------------------------------------------------------------------------------------
#definiowanie ustawien solwera:
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]		
	),
	#silnik intgracji rownan ruchu Newtona:
	NewtonIntegrator(damping=.4),	
	PyRunner(iterPeriod=1,command='letSave()'),
	PyRunner(command='letPlot()',iterPeriod=100),
        PyRunner(iterPeriod=100,command='print voxelPorosity(100,mincorner,O.bodies[3].state.pos)'),
	]

def letSave():
   if O.iter == maxIter-1:
     from yade import export
     export.text(file+'_'+str(number)+".yade")
     print "---------------------------"
     print "maxIter = ",maxIter
     print "mumber of particles = ",len(O.bodies)
     print "hight current = ", O.bodies[3].state.pos[1]-thickness/2.0
     print "porosity = ",utils.porosity(length*(O.bodies[3].state.pos[1]-thickness/2.0)*width) 

def letPlot():
    plot.addData(step=O.iter,e=utils.porosity(length*(O.bodies[3].state.pos[1]-thickness/2.0)*width) )

plot.plots={'step':('e')}
plot.plot()

from yade import qt; qt.Controller();qt.View()

#definiowanie warunkow kompresji:
O.engines = O.engines+[KinemCTDEngine(compSpeed=v,sigma_save=(),temoin_save=(),targetSigma=40000.0,LOG=False)]
O.dt = 0.4*PWaveTimeStep()

#???
triax.goal1 = triax.goal2 = triax.goal3 = 10000

import sys
while triax.porosity > targetPorosity:
  spFrAngle = 0.95*spFrAngle
  setContactFriction(radians(spFrAngle))
  print "\r Friction:",spFrAngle
  print "\r Porosity:",triax.porosity  
  sys.stdout.flush()	#???
  O.run(500,1)

from yade import plot, utils

print "\r --------"
print "\r Final friction:",spFrAngle
print "\r Final porosity:",triax.porosity

from yade import export
export.text(file+'.yade')

O.pause()

#--------------------------------------------------------------------------------------------

Regards,

Wojtek

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.