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 => Solved

Luc Scholtès confirmed that the question is solved:
Well, you are right, I made a mistake in my script when asking for the
Qin and Qout. It gives the correct result if, of course, you ask for the
right ids (...).

Nonethless, this is not what caused the problem initially and I am still
facing this problem with the solver. Unfortunately, you or Bruno cannot
obtain the same error with my scripts and that's a big mystery for me at
the moment... Since another person seems to face the same problem (so
this is not just a curse), we need to figure out if this is really a
problem with the solver or the scripts... I must say that I am confident
in saying that this is related to the solver in association with the
computer set up since I don't have any problem when I run my test with
useSolver=0.

I put here again the script that gives the following behaviors when
useSolver=3:

1) when the packing is not imported from a text file:
- works every time with yade installed with sudo apt-get install yade ( VTK file corrupted)
- works randomly with yade installed from  git sources if no other simulation is running at the same time
- never works with yade installed from  git sources if another simulation is running at the same time

2) when the packing is imported from a text file:
- works every time with yade installed with sudo apt-get install yade (VTK file corrupted)
- works randomly with yade installed from  git sources if no other simulation is running at the same time
- never works with yade installed from  git sources if another simulation is running at the same time

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/40.,returnSpherePack=True)
sp.toSimulation(material=sphereMat)

## if you uncomment, it imports a pre-existing packing
#PACKING='111_10k'
#O.bodies.append(ymport.text(PACKING+'.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

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)
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=0
        ,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(ids[0]) 
Qout = flow.getBoundaryFlux(ids[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.