Question #659557 on Yade changed:
https://answers.launchpad.net/yade/+question/659557
Status: Answered => Open
Luc Scholtès is still having a problem:
Here you go. It fails every time. Of course, you'll have to create your
own packing (cf PACK parameter in the first lines).
Also, I tried to simulate an injection and it does not look good at all
compared to what I used to obtain. Any chance you would have a MWE that
I could try on my side?
Thanks for your help.
Luc
from yade import pack, ymport
#### packing you want to test
PACK='111_5k'
intR=1.262
#### if you want to adjust the macroscopic permeability
pFactor=8.e-18
# for k=1e-20 m2: 8.e-18 for 111_5k
#### material definition
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))
#### preprocessing to get dimensions
O.bodies.append(ymport.text(PACK+'.spheres',scale=1.,shift=Vector3(0,0,0),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
R=0
Rmax=0
Rmin=1000
numSpheres=0.
for o in O.bodies:
if isinstance(o.shape,Sphere):
numSpheres+=1
R+=o.shape.radius
if o.shape.radius>Rmax:
Rmax=o.shape.radius
if o.shape.radius<Rmin:
Rmin=o.shape.radius
Rmean=R/numSpheres
#print "nb spheres=",numSpheres, " | mean Diameter=", 2*Rmean,
",X,Y,Z,", X,Y,Z
#### reset the scene now
O.reset()
#### simulation is built up now
mn,mx=Vector3(xinf,yinf,zinf),Vector3(xsup,ysup,zsup)
walls=aabbWalls([mn,mx],oversizeFactor=1.5,thickness=X/10.,material=wallMat)
wallIds=O.bodies.append(walls)
O.bodies.append(ymport.text(PACK+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
flow=DFNFlowEngine(
isActivated=1,
useSolver=3,
boundaryUseMaxMin = [0,0,0,0,0,0],
permeabilityFactor=pFactor,
viscosity=0.001
)
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(smoothJoint=True)]
)
,flow
,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
,NewtonIntegrator(damping=0.4)
]
## if you want to block the particles (better for permeability test)
for o in O.bodies:
o.state.blockedDOFs+='xyzXYZ'
## flow along X direction
print 'Flow along X!'
flow.bndCondIsPressure = [1,1,0,0,0,0]
flow.bndCondValue = [1,0,0,0,0,0]
O.step()
## getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(0)
Qout = flow.getBoundaryFlux(1)
# if Qout is the total discharge, we can compute k=Q*nu*Length/(Area*(Pout-Pin))
# if Qout is the flux, we can compute k=Q*nu*Length/(Pout-Pin) -> getFlux gives
total discharge -> Qout (m3/s)!
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() # if you want to see the result in Paraview
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
