yade-users team mailing list archive

[Othman Sh <question688195@xxxxxxxxxxxxxxxxxxxxx>]

Question #688195 on Yade changed:
https://answers.launchpad.net/yade/+question/688195

Description changed to:
Hello,

I'm trying to simulate 1-D flow through a cylindrical sphere packing
using the PFV model (FlowEngine) in Yade to mimic water flow through
porous concrete . So I want the spheres position to be fixed. I want to
impose a pressure on the top surface of the cylinder and get the flux at
the bottom surface. There should be no flux on the sides of the cylinder
(insulated, flow only in z direction). I have the following questions:

1- should the cylindrical sphere packing be inside a box? If I do the
simulation without the box I get a message like "Vh==NULL!! id=6
Point=0.08606 0.0994361 0.0201042 rad=0.0048993" which is repeated for
almost all body ids.

if the packing should be bounded by walls then I should have the packing
in a cylindrical wall, not in a box?

2- I don't know if my boundary conditions are correct for the case I'm trying to model. based on my understanding from the Yade documentation:
flow.bndCondIsPressure=[Xmin,Xmax,Ymin,Ymax,Zmin,Zmax] is that accurate?

so if I want my packing to have pressure on the top surface I should
have: flow.bndCondIsPressure=[0,0,0,0,0,1]. Correct?

In my code below (which I wrote it based on the Oedometer code by B.
Chareyre), the reason I have my boundary like this:
flow.bndCondIsPressure=[0,0,0,0,1,1] is that in the documentation it
says: "bndCondIsPressure: defines the type of boundary condition for
each side. True if pressure is imposed, False for no-flux". So if I want
water to pass the bottom surface (zmin) and to get the flux at zmin I
think I should have zmin turned to 1 in the  "flow.bndCondIsPressure
[]", correct?

I appreciate any help.

Thank you,
Othman

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


from yade import pack, ymport, plot

num_spheres=3000# number of spheres
young=1e6
radiuscyl=.05
heightcyl=.203
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.101,0.101,0.2035) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0049,rRelFuzz=0.002,seed=1)

					##	
#### cylinder extraction 										
pred=pack.inCylinder((.05,.05,0),(.05,.05,heightcyl),radiuscyl) 						
spFilter=filterSpherePack(pred,sp,returnSpherePack=True)				
spFilter.toSimulation(material='spheres')

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


flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,10]
flow.boundaryUseMaxMin=[0,0,0,0,1,1]
O.dt=0.1e-3
O.dynDt=False


O.run(1,1)
Qin = flow.getBoundaryFlux(5)# 5 is top wall 
Qout = flow.getBoundaryFlux(4)
Qxmin=flow.getBoundaryFlux(0)
Qxmax=flow.getBoundaryFlux(1)
Qymin=flow.getBoundaryFlux(2)
Qymax=flow.getBoundaryFlux(3)


NewtonIntegrator(damping=0)

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
## a function saving variables
def history():
	    plot.addData(Qin=Qin,Qout=Qout,t=O.time,p=flow.averagePressure(),Qxmin=Qxmin,Qxmax=Qxmax,Qymin=Qymin,Qymax=Qymax)
	plot.saveDataTxt('1.txt')


plot.plots={'t':('p')}
plot.plot()

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