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://answers.launchpad.net/yade/+question/688685
[2] wireframe: https://drive.google.com/open?id=1nx5dLhYYd0g0hF-T-j74gs4mnXv9pSam
[3] surface: https://drive.google.com/open?id=1Ftug3740yxHuuE0JJvsb71xafDoZJKUy
[4] wireframe: https://drive.google.com/open?id=13e2Pq9pZE3phjqBgqKI7MswkfsMxMzit
[5] surface: https://drive.google.com/open?id=1prgBJPlr5jJGrlhOnyj2q08Ot8s7TUzS
[6] https://drive.google.com/open?id=17KEDSWCspG8SJtd4B3U32YxMhD-wzjm0

------------------------------------------------

# -*- 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,allz,r=np.genfromtxt('30Porcyl.txt', unpack=True) #access the packed cyl file (txt)
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,mnz,mxx,mxy,mxz)
99998))

young=1e6

O.materials.append(FrictMat(young=young,poisson=0.2,density=1900,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz) # corners of the box that contain the packign. Obtained from aabbExtrema
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

############################# spheres #############################
spheres=O.bodies.append(ymport.text('/home/auser/DEM_codes/Flow/30Porcyl.txt'))

Height=max(utils.aabbDim())

#Fix all particles in their positions. No deformation
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"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 NewtonIntegrator(damping=0.2)
]

#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.permeabilityFactor=1
flow.viscosity=visc
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,dP]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e-6
O.dynDt=False

O.run(1,1)

################################################ blockcells ############################################################
########################################################################################################################

numPoints = 100
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints)
ys = np.linspace(0.95*mny,1.05*mxy,numPoints)
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)

cellsHit = []
#get the point at the center of the cylinder cross-section
cylcenterx=(mxx+mnx)/2
cylcentery=(mxy+mny)/2

for x,y,z in itertools.product(xs, ys, zs):
 cellId = flow.getCell(x,y,z) #getCell return cell id for a point in xyz
 if cellId in cellsHit: continue
 if np.sqrt((x-cylcenterx)**2+(y-cylcentery)**2) > radiuscyl: #this is to isolate the cells that are outside the cylinder. radiuscyl is the cylinder radius
  cellsHit.append(cellId)

for i in cellsHit:
# flow.blockCell(i,blockPressure=True) ##I tried this function but it didn't work and the cells in cellsHit were still counted in the flow model, so I used setCellPressure instead.
 flow.setCellPressure(i,0)
O.run(1,1)
qq=flow.averageVelocity()
Q=abs(qq[2])
Permeability=abs((density*g*Q*Height)/dP) #hydraulic conductivity in distance/time
print ("average velocity: ",Q, "Permeability m/s: ",Permeability,"Permeability in/hr: ", Permeability*141732)

########################################################################################################################
########################################################################################################################
print ("runtime = ", time.time()-tic)
#O.run()
#O.engines=O.engines+[PyRunner(iterPeriod=20,command='flow.saveVtk(withBoundaries=True)')]