permeability - FlowEngine
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:/
[2] https:/
[3] https:/
[4] This is a clip showing the vertical cross-section of the box https:/
[5] https:/
[6] https:/
#######
# -*- coding: utf-8 -*-
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=1e-3 #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,
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,
young=1e6
O.materials.
O.materials.
#create walls
mn,mx=Vector3(
walls=aabbWalls
wallIds=
yade.qt.View()
#######
spheres=
Height=
dP=P/Height #pressure gradient Pa/m
#Fix all particles in their positions. No deformation
for i in O.bodies:
i.state.
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
FlowEngine(
GlobalStiffnes
NewtonIntegrat
]
flow.dead=0
flow.useSolver=3
flow.permeabili
flow.viscosity=visc
flow.bndCondIsP
flow.bndCondVal
flow.boundaryUs
O.dt=1e-6
O.dynDt=False
flow.emulateAct
#######
numPoints = 100
xs = np.linspace(
ys = np.linspace(
zs = np.linspace(
cellsHit = [] #create array
cylcenterx=
cylcentery=
for x,y,z in itertools.
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(
cellsHit.
for i in cellsHit:
flow.blockCell
flow.setCellPr
O.run(1,1)
flow.meshUpdate
flow.defToleran
qq=flow.
Q=abs(qq[2]) #this is wrong, you need average velocity for the cylinder only not the whole system
Permeability=
print ("average velocity: ",Q, "Permeability m/s: ",Permeability,
#######
Z_height=mxz*0.95
cylcenterx=
cylcentery=
numPoints = 20
xr = np.linspace(
yr = np.linspace(
zr = np.linspace(
radiuscyl=
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.
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(
cellsHit2.
velocityVector = np.array(
velMag = np.linalg.
cellVol = flow.getCellVol
v = v + cellVol*
totalVolume += cellVol #total volume = total volume + cell vol
Cylinder_Vel = np.linalg.
Permeability2=
print("Cylinder_Vel : ",Cylinder_
print ("runtime = ", time.time()-tic)
#O.run()
#O.engines=
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Othman Sh
- Solved:
- Last query:
- Last reply: