yade-users team mailing list archive

[Othman Sh <question688685@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello All,

I wrote a code that is supposed to simulate a permeability test of a saturated sphere packing in a column. The spheres packing is rigid and no deformation/change in particle positions due to fluid flow should occur. However, I'm getting errors such as:
- Vh==NULL!!, 
- negative volume for an ordinary pore (temp warning, should still be safe)

Can anyone help me to make this code work? My goal is to get the average fluid velocity in the pores. The code is copied below. 

Thank you,

Othman

---------------------------
# -*- coding: utf-8 -*-

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


radiuscyl=.05
heightcyl=.2

dP=1e3 #Pa
visc=1e-3 #Pa.sec
density=1000 #kg/m3

O.materials.append(FrictMat(young = 5e10, poisson = 0.15,frictionAngle = atan(.2), density=1920))

############################ spheres #############################
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.3,.3,2),rMean=0.0083,rRelFuzz=0.1)	
#### cylinder extraction 										
pred=pack.inCylinder((.2,.2,0),(.2,.2,heightcyl),radiuscyl) 						
spFilter=filterSpherePack(pred,sp,Material=Material, returnSpherePack=True)				
spFilter.toSimulation()
############################ facets #############################

facets=geom.facetCylinder((.2,.2,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=150,wallMask=4)
cylinder=O.bodies.append(facets)
yade.qt.View()
Height=max(utils.aabbDim())

#Fix all particles in their positions. No deformation
for i in O.bodies:
	i.state.blockedDOFs='xyzXYZ'

############################ Engines #############################
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"),#introduced as a dead engine for the moment, see 2nd section
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	NewtonIntegrator(damping=0.2)
	
]
O.dt=1e-5
O.run(1,1)

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.useSolver=3
flow.viscosity=visc
flow.bndCondIsPressure=[0,0,0,0,1,1] #[xmin,xmax,ymin,ymax,zmin,zmax]
flow.bndCondValue=[0,0,0,0,0,dP]
flow.boundaryUseMaxMin=[0,0,0,0,1,1] #if you want wall to be boundary, set this 0. If you don't have a wall, make it 1 so that spheres will be boundary


#Permeability calculations
Q=flow.averageVelocity()
Permeability=visc*Q*Height/dP
print ("average velocity: ",Q, "Permeability: ",Permeability)

O.run()

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