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