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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
