# PFV model

Hello,

I'm trying to obtain the permeability of spheres packed in a cylindrical shape using the PFV model in Yade. I previously posted a question about this problem [1]. The PVF only worked when I have my cylinder inside a box. The problem is that the box corners affects the flow results. Robert Caulk provided a solution to get the permeability in specific volume inside the cylinder [1, #20]. I tried that but the flow velocity results were high and unrealistic. I tried something that was suggested by Bruno Chareyre [1, #19] which is to block/deactivate the cells that are outside the cylinder core so that they're not included in the flow engine. I did that and I used Paraview to see the results and the cells around the cylinder were blocked, but the applied pressure seems to be only affecting the top surface [2, 3]. However, after the first ~1000 run, the cells were active again and the velocity suddenly jumps from 0.087 to 2269, you can see the Paraview visualization [4, 5]. Any idea why this happens and is there something wrong in my code (copied below)?

Thank you,

Othman

P.S. to run the code, the sphere packing txt file can be found in [6].

[1] https:/

[2] wireframe: https:/

[3] surface: https:/

[4] wireframe: https:/

[5] surface: 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

heightcyl=.25

dP=4.347e3 #Pa

visc=1e-3 #Pa.sec - taken from Catalano, Chareyre, Bathelemy (2014)

density=1000 #kg/m3

g=9.81 #m/s2

allx,ally,

mnx=min(allx)*0.99 #to get the xyz of the lowest corner, multiply by factor to reduce it to let the walls surround all the spheres.

mny=min(ally)*0.99

mnz=min(allz)*0.99

mxx=max(allx)*1.01 #to 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,

99998))

young=1e6

O.materials.

O.materials.

mn,mx=Vector3(

walls=aabbWalls

wallIds=

#######

spheres=

Height=

#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

]

#Activate flow engine and set boundary conditions in order to get permeability

flow.dead=0

flow.useSolver=3 #CHOLMOD: sparse direct solver for matrices. Its a library that is part of SuiteSparse linear algebra package.

flow.permeabili

flow.viscosity=visc

flow.bndCondIsP

flow.bndCondVal

flow.boundaryUs

O.dt=1e-6

O.dynDt=False

O.run(1,1)

#######

#######

numPoints = 100

xs = np.linspace(

ys = np.linspace(

zs = np.linspace(

cellsHit = []

#get the point at the center of the cylinder cross-section

cylcenterx=

cylcentery=

for x,y,z in itertools.

cellId = flow.getCell(x,y,z) #getCell return cell id for a point in xyz

if cellId in cellsHit: continue

if np.sqrt(

cellsHit.

for i in cellsHit:

# flow.blockCell(

flow.setCellPr

O.run(1,1)

qq=flow.

Q=abs(qq[2])

Permeability=

print ("average velocity: ",Q, "Permeability m/s: ",Permeability,

#######

#######

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

- Last query:
- 2020-03-05

- Last reply:
- 2020-03-05

Robert Caulk (rcaulk) said : | #1 |

>The problem is that the box corners affects the flow results.

Please elaborate on how the flow results are affected.

>However, after the first ~1000 run, the cells were active again and the velocity suddenly jumps from 0.087 to 2269

>I tried this function but it didn't work and the cells in cellsHit were still counted in the flow model

Try blocking the cells before the mesh is created (before a O.run(1,1) where FlowEngine is active).

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

Hi Robert,

>Please elaborate on how the flow results are affected.

The cylinder is inside a box, so the box corners are empty (no spheres) so the flow velocity will be very high. In the code that I posted, if I comment everything after "blockcells" section, the averageVelocity() in z direction = 2269.56 which is super high (it should be less than 1). Also the visualization in Paraview shows the difference between modeling the flow only in the cylinder [1] and when the whole box is included in the model [2].

>Try blocking the cells before the mesh is created (before a O.run(1,1) where FlowEngine is active).

I did that but I get this error: Triangulation does not exist. Sorry.

Its weird that the code works for 1000 run then suddenly the cells in cellsHit are not blocked anymore.

[1] https:/

[2] https:/

Robert Caulk (rcaulk) said : | #3 |

Try setting

flow.meshUpdate

flow.defToleran

to avoid remeshing.

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

That worked really well!

Thanks, Robert. You are awesome!. Problem solved