filtration - PeriodicFlowEngine

Dear all,
I would like to simulate the transport of small particles by downward flow like was done in (Sari et al. 2011). My packing has periodic lateral boundaries; a bottom plane and another top plane are introduced to the simulation.
It seems that the PeriodicFlowEngine is not working correctly (the fine particle doesn't move), I am wondering if the setting for this engine were chosen properly (flow.bndCondIsPressure,flow.bndCondValue, flow.gradP), or if the fine particle at the top of the packing may influence the triangulation in this case.
Does anyone know what I did wrong? And how can I fix that ?
Thank you!
Feda

Sari, H., Chareyre, B., Catalano, E., Philippe, P., & Vincens, E. (2011). Investigation of internal erosion processes using a coupled dem-fluid method. In Particles 2011 II International Conference on Particle-Based Methods, E. Oate and DRJ Owen (Eds), Barcelona (pp. 1-11).

from yade import pack,qt
import math
O.periodic=True

nbofspheres=1500
width=0.04
depth=0.04
height=0.1
radius=0.002
radiusfine=0.0005
rRelFuzz=0.01

O.cell.hSize=Matrix3(width, 0, 0,
            0 ,depth , 0,
              0, 0, 3*height)

O.materials.append(ViscElMat(kn=1e7,ks=1e7,en=0.1,et=0.1,frictionAngle=0,density=2530,label='spheres'))
O.materials.append(ViscElMat(kn=1e7,ks=1e7,en=0.1,et=0.1,frictionAngle=0,density=2530,label='walls'))

bottomWall = utils.box(center=(width/2,depth/2,height-0.01), extents=(width*1000,depth*1000,0) ,fixed=True,wire=False,material='walls',color=(0,0,1))
topWall = utils.box(center=(width/2,depth/2,2*height+0.01), extents=(width*1000,depth*1000,0) ,fixed=True,wire=False,material='walls',color=(0,0,1))
O.bodies.append([bottomWall,topWall])

sp=pack.SpherePack()
sp.makeCloud((0,0,height),(width,depth,2*height),radius,rRelFuzz,nbofspheres,periodic=True,porosity=0.3,seed=1)
O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,1),material='spheres') for s in sp])
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.state.blockedDOFs='xyzXYZ'
        b.state.vel=(0,0,0)
        b.state.angVel=(0,0,0)

O.bodies.append(utils.sphere((random.uniform(0,0.04),random.uniform(0,0.04),2*height+0.005),radiusfine,highlight=True,material='spheres',color=(1,0,0)))
sp.fromSimulation()

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],allowBiggerThanPeriod=True),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
 [Law2_ScGeom_ViscElPhys_Basic()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=0.8*PWaveTimeStep(),viscEl=True),
    PeriodicFlowEngine(dead=1,label='flow'),
    NewtonIntegrator(damping=0,gravity=(0,0,0),label='newton')
    #PyRunner(command='checkPositions()',realPeriod=10),
]

flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=1
flow.wallIds=[-1,-1,-1,-1,1,0]
flow.bndCondIsPressure=[0,0,0,0,1,0]
flow.bndCondValue=[0,0,0,0,-50000,0]
flow.boundaryUseMaxMin=[0,0,0,0,1,1]
#flow.gradP=(0,0,-50000)

O.run()
qt.View()