yade-users team mailing list archive

[Lingran Zhang <question686131@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello,

One question for FlowEngine and DFNFlowEngine: why the flow.nCells() != "cell number" recorded by flow.saveVtk()?

I would like to record the Ids of initial fractured cells and then impose different pressures within them, but I find that the ids (and vertices) are not consistent given by these two ways.

Below is an example.
Thank you in advance for your help.
Best regards,
Lingran


from yade import pack, ymport

### material
def sphereMat(): return JCFpmMat(type=1,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(18))
def wallMat(): return JCFpmMat(type=0,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(0))

### bodies
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
O.bodies.append(aabbWalls([mn,mx],oversizeFactor=1.5,thickness=0.1,material=wallMat))
sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(1,1.,1.)),spheresInCell=100,radius=1/10.,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat)

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

### engines
flow=FlowEngine(
        isActivated=1,
        ### choose solver to use (0: Gauss Seidel, 1: Taucs, 2: Pardiso, 3: CHOLMOD)
        useSolver=3, # 3 should be used by default but does not work everytime...
         boundaryUseMaxMin = [0,0,0,0,0,0],
        bndCondIsPressure = [1,1,0,0,0,0],
        bndCondValue = [1,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        #debug=1
)

intR=1.1
O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')])
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8)
 ,NewtonIntegrator()
]

### simulation starts here

## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

O.run(1,True)
## getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(0)
Qout = flow.getBoundaryFlux(1)


permeability = abs(Qout)*flow.viscosity*X/(Y*Z) # !!! if Pout=1, Pin=0
permeability2 = flow.averageVelocity()*flow.viscosity*X # !!! if Pout=1, Pin=0
conductivity = permeability*1000*9.82/flow.viscosity # K=rho*g*k/nu
print "Qin=",Qin," Qout=",Qout," ARE THEY EQUAL? IF NOT-> NO FLOW!"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # to see the result in Paraview


for i in range(flow.nCells()):
   print i, flow.getVertices(i)

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