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:
- OK for the boundaryUseMaxMin=[1,1,1,1,1,1] . I misunderstood the doc.

- For the calibration of the coordination number, yes, we could do that
with a dedicated python script but I must say that I am not a big fan of
the randomDensePack() function, even though my way to generate dense
packing is probably not better.

- Regarding the problem with the solver, I guess I'll need to figure it
out by myself... Not an easy task...

- Regarding the walls appended after the packing, no error message shows
up really (even with debug=1). It is just that the computation is wrong
(Qin and Qout are not equal) and that the VTK file cannot be open with
paraview (error reading ascii cell data). Could this be related to the
tesselation?

Here is the script:

from yade import pack, ymport

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

#### bodies
X=Y=Z=1

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(X,Y,Z)),spheresInCell=1000,radius=X/20.,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

Rmean=0
nbSph=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
Rmean=Rmean/nbSph

print "nb spheres=",nbSph, " | mean Diameter=", 2*Rmean, ",X,Y,Z=",
X,Y,Z

mn,mx=Vector3(xinf,yinf,zinf),Vector3(xsup,ysup,zsup)
#mn,mx=Vector3(0,0,0),Vector3(X,Y,Z)
walls=utils.aabbWalls(extrema=[mn,mx],oversizeFactor=1.5,thickness=min(X,Y,Z)/100.,material=sphereMat)
ids=O.bodies.append(walls)
print 'walls ids=',ids

#### engines
flow=FlowEngine(
        isActivated=1,
        useSolver=3,
        wallIds=ids,
        boundaryUseMaxMin=[0,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        debug=1
)

O.engines=[
	ForceResetter()
	,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()])
	,InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
	)
	,flow
	,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
	,NewtonIntegrator(damping=0.)
]

#### simulation
## 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?"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # if you want to see the result in Paraview

If you want to compare with the script where the walls are appended
before the packing (which gives a reasonable result with Qin=Qout and a
vtk file that can be open in Paraview), here it is:

from yade import pack, ymport

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

#### bodies
X=Y=Z=1

mn,mx=Vector3(0,0,0),Vector3(X,Y,Z)
walls=utils.aabbWalls(extrema=[mn,mx],oversizeFactor=1.5,thickness=min(X,Y,Z)/100.,material=sphereMat)
ids=O.bodies.append(walls)
print 'wall ids=',ids

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(X,Y,Z)),spheresInCell=1000,radius=X/20.,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

Rmean=0
nbSph=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
Rmean=Rmean/nbSph

print "nb spheres=",nbSph, " | mean Diameter=", 2*Rmean, ",X,Y,Z=",
X,Y,Z

#### engines
flow=FlowEngine(
        isActivated=1,
        useSolver=3,
        wallIds=ids,
        boundaryUseMaxMin=[0,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        debug=1
)

O.engines=[
	ForceResetter()
	,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()])
	,InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
	)
	,flow
	,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep())
	,NewtonIntegrator(damping=0.)
]

#### simulation
## 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?"
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.