yade-users team mailing list archive

[Wojciech Sobieski <question255937@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Wojciech Sobieski is still having a problem:
Hi Bruno!

1. I use a slightly different way. First I determine the target porosity
and the number of particles, and next I calculate the theoretical volume
of the bed. It turned out that the final volume is bigger as the initial
volume of the cloud (given by min and max corners). After several
attempts I found the optimal size of the initial volume. The defect is
that the average radius is in every case is a little different. For this
reason I want to know how to calculate the average radius in YADE.

2. I tried likewise to resolve my task in this way. The problem is that
the calculation time is very long and in fact I don't understand fully
the python script. I tried to connect the trick with precisely defined
target porosity and the example with triax test. I afraid that some
parts may be without sense. For example I still don't know what is the
triax.goal.

If you have few minutes - check my script below and give me some remarks
how to improve it, please. If not, I will try do it alone.

Regards,

Wojtek.

-------
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 = 1000

#predkosc kompresji (scianki, ktora sie rusza):
v = 1.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 

#--------------------------------------------------------------------------------------------
#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) 

#kasowanie danych do wizualizacji:
plot.resetData()

#--------------------------------------------------------------------------------------------
#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()]		
	),
	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=.4*PWaveTimeStep()

triax.goal1 = triax.goal2 = triax.goal3 = 10000

import sys #this is only for the flush() below
while triax.porosity > targetPorosity:
  compFricDegree = 0.95*compFricDegree
  setContactFriction(radians(compFricDegree))
  print "\r Friction:",compFricDegree
  print "\r Porosity:",triax.porosity  
  sys.stdout.flush()
  # while we run steps, triax will tend to grow particles as the packing
  # keeps shrinking as a consequence of decreasing friction. Consequently
  # porosity will decrease
  O.run(500,1)

from yade import plot, utils

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

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

O.pause()

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