PFV model

Asked by Othman Sh on 2020-03-03

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)')]

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://drive.google.com/file/d/1Ftug3740yxHuuE0JJvsb71xafDoZJKUy
[2] https://drive.google.com/open?id=1prgBJPlr5jJGrlhOnyj2q08Ot8s7TUzS

Robert Caulk (rcaulk) said : #3

Try setting

flow.meshUpdateInterval=-1
flow.defTolerance=-1

to avoid remeshing.

Othman Sh (othman-sh) said : #4

That worked really well!

Thanks, Robert. You are awesome!. Problem solved