yade-users team mailing list archive

[Luc Scholtès <question659557@xxxxxxxxxxxxxxxxxxxxx>]

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.