yadeusers team mailing list archive

yadeusers team

Mailing list archive

Message #22512
[Question #689330]: permeability  FlowEngine
New question #689330 on Yade:
https://answers.launchpad.net/yade/+question/689330
Hello,
I previously asked about PFV FlowEngine and posted several questions about it [1,2]. I'm trying to get the permeability of a cylindrical sphere packing. The PVF engine only worked when I have my cylinder inside a box [1, #13,14,16]. Since the cylinder is inside a box, PFV will create cells in all the box (not just the cylinder) and the box corners are empty (no spheres) so the flow velocity will be very high. I am trying to block the cells that are outside the cylinder, it seems that this is not working. Using Paraview, I can see that the pressure is 0 outside the cylinder but the velocity is not. I am using blockCell() but I'm not sure I am using it proberly. Using this should exclude the cells from the PFV model so the velocity in those excluded cells should be 0, right? That is what I understood from Yade FlowEngine documentation. I exported the results in VTK with borders and in [3] and [4], you can see the velocity values are not zero in the corners.
Just to clarify, the reason I want to have the cylinder "sealed" from the sides is to mimic the lab experiment that I am doing in which I have my cylinder wrapped with insulation material and the flow is onle from top surface to the bottom (no lateral flow from the sides) [5].
Thanks,
Othman
P.S. to run the code, the sphere packing txt file can be found in [6]. This makes running the code much faster than packing spheres every run.
[1] https://answers.launchpad.net/yade/+question/688685
[2] https://answers.launchpad.net/yade/+question/689117
[3] https://drive.google.com/open?id=1_SRG2bWYZuLaRKU5MzuV3UzJ5otAtbj
[4] This is a clip showing the vertical crosssection of the box https://drive.google.com/open?id=1IEagKzpZgRCRVksGj1SrpK2v9cpCoM_6
[5] https://drive.google.com/open?id=1plzJ6SVpgghihhb7N3lLAPqrOR0Ra0j
[6] https://drive.google.com/open?id=17KEDSWCspG8SJtd4B3U32YxMhDwzjm0
##################################################
# * coding: utf8 *
from builtins import range
from yade import pack, ymport
import numpy as np
import time
tic=time.time()
radiuscyl=.05 #cylinder radius
P=4500 # pressure at top surface in Pa
visc=1e3 #Pa.sec  taken from Catalano, Chareyre, Bathelemy (2014)
density=1000 #kg/m3
g=9.81 #m/s2
# to access the sphere packing and get the minimum and maximum corner points
allx,ally,allz,r=np.genfromtxt('30Porcyl.txt', unpack=True) #access the packed cyl file (txt)
mnx=min(allx)*0.99 #same as aabbExtrema, get the xyz of the lowest corner, multiply by factor to reduce it to let the wall surround all the spheres.
mny=min(ally)*0.99
mnz=min(allz)*0.99
mxx=max(allx)*1.01 #same as aabbExtrema, get the xyz of the top corner, multiply by factor to increase it to let the wall surround all the spheres.
mxy=max(ally)*1.01
mxz=max(allz)*1.01
print ('mn xyz, mx xyz', mnx,mny,mnz,mxx,mxy,mxz)
young=1e6
O.materials.append(FrictMat(young=young,poisson=0.2,density=1900,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
#create walls
mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
yade.qt.View()
############################# spheres #############################
spheres=O.bodies.append(ymport.text('/home/auser/DEM_codes/Flow/30Porcyl.txt'))
Height=max(utils.aabbDim())
dP=P/Height #pressure gradient Pa/m
#Fix all particles in their positions. No deformation
for i in O.bodies:
i.state.blockedDOFs='xyzXYZ'
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()],label="iloop"
),
FlowEngine(dead=1,label="flow"),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
NewtonIntegrator(damping=0.2)
]
flow.dead=0
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=visc
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,P]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e6
O.dynDt=False
flow.emulateAction()
############################ blockcells #######################
numPoints = 100
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints) #np.linspace is to divide the distance between 2 points
ys = np.linspace(0.95*mny,1.05*mxy,numPoints) #if the factor for mx change to less than 1, code will not work properly.
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)
cellsHit = [] #create array
cylcenterx=(mxx+mnx)/2
cylcentery=(mxy+mny)/2
for x,y,z in itertools.product(xs, ys, zs):
cellId = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
if cellId in cellsHit: continue #this prevents counting a cell multiple times
if np.sqrt((xcylcenterx)**2+(ycylcentery)**2) > radiuscyl:
cellsHit.append(cellId)
for i in cellsHit:
flow.blockCell(i,blockPressure=True)
flow.setCellPressure(i,0)
O.run(1,1)
flow.meshUpdateInterval=1 #these two lines to prevent remeshing after 1000 run which unblock all cells in cellsHit
flow.defTolerance=1
qq=flow.averageVelocity()
Q=abs(qq[2]) #this is wrong, you need average velocity for the cylinder only not the whole system
Permeability=abs((density*g*Q*Height)/dP)
print ("average velocity: ",Q, "Permeability m/s: ",Permeability,"Permeability in/hr: ", Permeability*141732)
######################## Calculate average velocity at the cylinder ##########################
Z_height=mxz*0.95
cylcenterx=(mxx+mnx)/2
cylcentery=(mxy+mny)/2
numPoints = 20
xr = np.linspace(0.95*mnx,1.05*mxx,numPoints) #np.linspace is to divide the distance between 2 points
yr = np.linspace(0.95*mny,1.05*mxy,numPoints)
zr = np.linspace(0.95*mnz,Z_height,numPoints)
radiuscyl=radiuscyl*0.98
totalVolume = 0 #zero a variable
v = np.array([0,0,0]) #what is the difference between this and cellhit array
##get the cellIDs for a given portion of the cylinder
cellsHit2=[]
for x,y,z in itertools.product(xr, yr, zr):
cellId2 = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
if cellId in cellsHit2: continue #what does this mean? CellsHit is empty before this line
if np.sqrt((xcylcenterx)**2+(ycylcentery)**2) < radiuscyl:
cellsHit2.append(cellId2)
velocityVector = np.array(flow.getCellVelocity((x,y,z))) #getCellVelocity is a tool created by Robert Caulk. It will return the vel vectors for the cell xyz
velMag = np.linalg.norm(velocityVector) #linalg return the magnitude of the velocity vector
cellVol = flow.getCellVolume((x,y,z,)) #getCellVolume is a tool created by Robert Caulk
v = v + cellVol*velocityVector
totalVolume += cellVol #total volume = total volume + cell vol
Cylinder_Vel = np.linalg.norm(v)/totalVolume
Permeability2=abs((density*g*Cylinder_Vel*Height)/dP)
print("Cylinder_Vel : ",Cylinder_Vel,"permeability 2 (in/hr): ",Permeability2*141732)
print ("runtime = ", time.time()tic)
#O.run()
#O.engines=O.engines+[PyRunner(iterPeriod=20,command='flow.saveVtk(withBoundaries=True)')]

You received this question notification because your team yadeusers is
an answer contact for Yade.