permeability model

Hello,

I'm trying to simulate 1-D flow through a cylindrical sphere packing using the PFV model (FlowEngine) in Yade to mimic water flow through porous concrete . So I want the spheres position to be fixed. I want to impose a pressure on the top surface of the cylinder and get the flux at the bottom surface. There should be no flux on the sides of the cylinder (insulated, flow only in z direction). I have the following questions:

1- should the cylindrical sphere packing be inside a box? If I do the simulation without the box I get a message like "Vh==NULL!! id=6 Point=0.08606 0.0994361 0.0201042 rad=0.0048993" which is repeated for almost all body ids.

if the packing should be bounded by walls then I should have the packing in a cylindrical wall, not in a box?

2- I don't know if my boundary conditions are correct for the case I'm trying to model. based on my understanding from the Yade documentation:
flow.bndCondIsPressure=[Xmin,Xmax,Ymin,Ymax,Zmin,Zmax] is that accurate?

so if I want my packing to have pressure on the top surface I should have: flow.bndCondIsPressure=[0,0,0,0,0,1]. Correct?

In my code below (which I wrote it based on the Oedometer code by B. Chareyre), the reason I have my boundary like this: flow.bndCondIsPressure=[0,0,0,0,1,1] is that in the documentation it says: "bndCondIsPressure: defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux". So if I want water to pass the bottom surface (zmin) and to get the flux at zmin I think I should have zmin turned to 1 in the "flow.bndCondIsPressure []", correct?

I appreciate any help.

Thank you,
Othman

-------------
# -*- coding: utf-8 -*-

from yade import pack, ymport, plot

num_spheres=3000# number of spheres
young=1e6
radiuscyl=.05
heightcyl=.203
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.101,0.101,0.2035) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0049,rRelFuzz=0.002,seed=1)

     ##
#### cylinder extraction
pred=pack.inCylinder((.05,.05,0),(.05,.05,heightcyl),radiuscyl)
spFilter=filterSpherePack(pred,sp,returnSpherePack=True)
spFilter.toSimulation(material='spheres')

yade.qt.View()
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"),#introduced as a dead engine for the moment, see 2nd section
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 NewtonIntegrator(damping=0.2)

]

flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,10]
flow.boundaryUseMaxMin=[0,0,0,0,1,1]
O.dt=0.1e-3
O.dynDt=False

O.run(1,1)
Qin = flow.getBoundaryFlux(5)# 5 is top wall
Qout = flow.getBoundaryFlux(4)
Qxmin=flow.getBoundaryFlux(0)
Qxmax=flow.getBoundaryFlux(1)
Qymin=flow.getBoundaryFlux(2)
Qymax=flow.getBoundaryFlux(3)

NewtonIntegrator(damping=0)

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
## a function saving variables
def history():
     plot.addData(Qin=Qin,Qout=Qout,t=O.time,p=flow.averagePressure(),Qxmin=Qxmin,Qxmax=Qxmax,Qymin=Qymin,Qymax=Qymax)
 plot.saveDataTxt('1.txt')

plot.plots={'t':('p')}
plot.plot()