# permeability model

Asked by Othman Sh on 2020-01-22

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
heightcyl=.203
compFricDegree = 3 # initial contact friction during the confining phase
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=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
spFilter=filterSpherePack(pred,sp,returnSpherePack=True)
spFilter.toSimulation(material='spheres')

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"
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
NewtonIntegrator(damping=0.2)

]

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.saveDataTxt('1.txt')

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

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Othman Sh
Solved:
2020-01-29
Last query:
2020-01-29
2020-01-29
 Othman Sh (othman-sh) said on 2020-01-25: #1

If someone can answer this question please do, I really appreciate any help.
Thanks

 Othman Sh (othman-sh) said on 2020-01-28: #2

Just following up on this question, any reply?
Thank you

 Robert Caulk (rcaulk) said on 2020-01-29: #3

Hello,

>>"Vh==NULL!! id=6 Point=0.08606 0.0994361 0.0201042 rad=0.0048993"

Decrease the time step.

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

That really depends on what you are trying to simulate. From an algorithmic perspective, it does not need walls (see boundaryUseMaxMin comment below).

>>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?

Yes.

>>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?

Yes.

>>flow.boundaryUseMaxMin=[0,0,0,0,1,1]

This command tells yade to use walls or max min coordinates to create the bounding spheres [2]. Right now, you are telling yade to use the wall locations for x y and use max min for z. Instead, use flow.boundaryUseMaxMin=[0,0,0,0,0,0] since you have 6 walls.

Cheers,

Robert

 Othman Sh (othman-sh) said on 2020-01-29: #4

Thanks Robert. That was really helpful and solves my problem.

Regards,
Othman