yade-users team mailing list archive

[Othman Sh <question689330@xxxxxxxxxxxxxxxxxxxxx>]

New question #689330 on Yade:
https://answers.launchpad.net/yade/+question/689330

Hello,

I previously asked about PFV FlowEngine and posted several questions about it [1,2]. I'm trying to get the permeability of a cylindrical sphere packing. The PVF engine only worked when I have my cylinder inside a box [1, #13,14,16]. Since the cylinder is inside a box, PFV will create cells in all the box (not just the cylinder) and the box corners are empty (no spheres) so the flow velocity will be very high. I am trying to block the cells that are outside the cylinder, it seems that this is not working. Using Paraview, I can see that the pressure is 0 outside the cylinder but the velocity is not. I am using blockCell() but I'm not sure I am using it proberly. Using this should exclude the cells from the PFV model so the velocity in those excluded cells should be 0, right? That is what I understood from Yade FlowEngine documentation. I exported the results in VTK with borders and in [3] and [4], you can see the velocity values are not zero in the corners. 

Just to clarify, the reason I want to have the cylinder "sealed" from the sides is to mimic the lab experiment that I am doing in which I have my cylinder wrapped with insulation material and the flow is onle from top surface to the bottom (no lateral flow from the sides) [5]. 

Thanks,
Othman 

P.S. to run the code, the sphere packing txt file can be found in [6]. This makes running the code much faster than packing spheres every run.


[1] https://answers.launchpad.net/yade/+question/688685
[2] https://answers.launchpad.net/yade/+question/689117
[3] https://drive.google.com/open?id=1_SRG2bWYZuLaRKU5Mzu-V3UzJ5otAtbj
[4] This is a clip showing the vertical cross-section of the box https://drive.google.com/open?id=1IEagKzpZgRCRVksGj1SrpK2v9cpCoM_6
[5] https://drive.google.com/open?id=1-plzJ6SVpgghihhb7N3lLAPqrOR0Ra0j
[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 #cylinder radius
P=4500 # pressure at top surface in Pa
visc=1e-3 #Pa.sec - taken from Catalano, Chareyre, Bathelemy (2014)
density=1000 #kg/m3
g=9.81 #m/s2

# to access the sphere packing and get the minimum and maximum corner points

allx,ally,allz,r=np.genfromtxt('30Porcyl.txt', unpack=True) #access the packed cyl file (txt)
mnx=min(allx)*0.99 #same as aabbExtrema, get the xyz of the lowest corner, multiply by factor to reduce it to let the wall surround all the spheres.
mny=min(ally)*0.99
mnz=min(allz)*0.99
mxx=max(allx)*1.01 #same as aabbExtrema, 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)


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'))
#create walls
mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz) 
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

yade.qt.View()

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

Height=max(utils.aabbDim())
dP=P/Height #pressure gradient Pa/m

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


flow.dead=0
flow.useSolver=3  
flow.permeabilityFactor=1
flow.viscosity=visc
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,P]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e-6
O.dynDt=False
flow.emulateAction()


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

numPoints = 100
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints) #np.linspace is to divide the distance between 2 points
ys = np.linspace(0.95*mny,1.05*mxy,numPoints) #if the factor for mx change to less than 1, code will not work properly. 
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)

cellsHit = [] #create array
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 position xyz
	if cellId in cellsHit: continue #this prevents counting a cell multiple times
	if np.sqrt((x-cylcenterx)**2+(y-cylcentery)**2) > radiuscyl:
		cellsHit.append(cellId)
		
for i in cellsHit: 
	flow.blockCell(i,blockPressure=True)
	flow.setCellPressure(i,0)
O.run(1,1)

flow.meshUpdateInterval=-1 #these two lines to prevent remeshing after 1000 run which unblock all cells in cellsHit
flow.defTolerance=-1

qq=flow.averageVelocity()
Q=abs(qq[2]) #this is wrong, you need average velocity for the cylinder only not the whole system
Permeability=abs((density*g*Q*Height)/dP)
print ("average velocity: ",Q, "Permeability m/s: ",Permeability,"Permeability in/hr: ", Permeability*141732)

######################## Calculate average velocity at the cylinder ##########################
Z_height=mxz*0.95 
cylcenterx=(mxx+mnx)/2
cylcentery=(mxy+mny)/2
numPoints = 20
xr = np.linspace(0.95*mnx,1.05*mxx,numPoints) #np.linspace is to divide the distance between 2 points
yr = np.linspace(0.95*mny,1.05*mxy,numPoints)
zr = np.linspace(0.95*mnz,Z_height,numPoints)
radiuscyl=radiuscyl*0.98 
totalVolume = 0 #zero a variable
v = np.array([0,0,0]) #what is the difference between this and cellhit array

##get the cellIDs for a given portion of the cylinder
cellsHit2=[]
for x,y,z in itertools.product(xr, yr, zr):
	cellId2 = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
	if cellId in cellsHit2: continue #what does this mean? CellsHit is empty before this line
	if np.sqrt((x-cylcenterx)**2+(y-cylcentery)**2) < radiuscyl:
	 cellsHit2.append(cellId2)
	 velocityVector = np.array(flow.getCellVelocity((x,y,z))) #getCellVelocity is a tool created by Robert Caulk. It will return the vel vectors for the cell xyz
	 velMag = np.linalg.norm(velocityVector) #linalg return the magnitude of the velocity vector
	 cellVol = flow.getCellVolume((x,y,z,)) #getCellVolume is a tool created by Robert Caulk
	 v = v + cellVol*velocityVector
	 totalVolume += cellVol #total volume = total volume + cell vol

Cylinder_Vel = np.linalg.norm(v)/totalVolume

Permeability2=abs((density*g*Cylinder_Vel*Height)/dP)
print("Cylinder_Vel : ",Cylinder_Vel,"permeability 2 (in/hr): ",Permeability2*141732)

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


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.