yade-users team mailing list archive

[Othman Sh <question693297@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi All,

I am trying to obtain the volume of individual pores in my sphere packing. I found this function in Yade documentation: 

volume( [ (int)id=0 ] ) → float
Returns the volume of Voronoi’s cell of a sphere

I have few questions here:

1- Is this the correct approach to get the pore volume for each individual cell ID?
2- Does this function gives the volume of the whole tetrahedron element (which includes part of the spheres surrounding the pore as in Fig.3a in [1]) or just the pore (i.e. fluid phase)? 
3- In my code below, this function is not working and I get this error: Segmentation fault (core dumped)
I have tried using this function in the terminal and I found that it sometimes give a value for some cell IDs and most of the time it will crush and this error occurs. My code is copied below. The packing txt file link is here [2]

I will appreciate your help in answering my questions. 

Thank you,
Othman


[1] https://onlinelibrary.wiley.com/doi/full/10.1002/nag.2198
[2] https://drive.google.com/file/d/1JqFe3W1cVFgKEa9lqU-MsCBxY7RDgY89/view?usp=sharing
---------------------------------------------

from yade import plot,pack, export, ymport
import numpy as np


P=4347 #Pa
visc=1e-3 #Pa.sec 
density=1000 #kg/m3
g=9.81 #m/s2
radiuscyl=.05
########## create walls ##########
allx,ally,allz,r=np.genfromtxt('S30.txt', unpack=True) #access the packed cyl file (txt)
mnx=min(allx)*0.999 
mny=min(ally)*0.999
mnz=min(allz)*0.999
mxx=max(allx)*1.001 
mxy=max(ally)*1.001
mxz=max(allz)*1.001

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0)
wallIds=O.bodies.append(walls)

########## spheres ##########
spheres=O.bodies.append(ymport.text('/home/auser/DEM_codes/Clogging/S30.txt')) #please change this when you download the packing file
yade.qt.View()
Height=max(utils.aabbDim())
dP=P/Height #Pa/m

for i in spheres:
	body= O.bodies[i]
	body.state.blockedDOFs='xyzXYZ'
print ('porosity = ', utils.porosity())

Height=max(utils.aabbDim())

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

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.useSolver=3 
flow.permeabilityFactor=1
flow.viscosity=visc
flow.bndCondIsPressure=[1,1,1,1,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()

#########################################################################################################################
################################################# Get pores volume#####################################
##########################################################################################################################
cellsHit = []
numPoints = 20
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)

v = np.array([0,0,0]) 


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
	cellsHit.append(cellId)
	VV=flow.volume(cellId)



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