# 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:
- 2020-03-19

- Last query:
- 2020-03-19

- Last reply:

Othman Sh (othman-sh) said : | #1 |

My problem was solved by changing the boundary condition from bndCondIsPressu