yade-users team mailing list archive

[Liu Changdong <question693485@xxxxxxxxxxxxxxxxxxxxx>]

New question #693485 on Yade:
https://answers.launchpad.net/yade/+question/693485

Hello everyone:
I am a beginner in PFV, i have a bit of problem with water pressure in PFV and DFNFlowEngine. I am trying to inject water in a  cubic rock matrix with flow.imposeFlux((0.5,0.5,0.5),-0.01). Although the program did not run incorrectly, I found that there was a large negative water pressure (-7.2e+7pa) in Paraview. So my question is : 
1. Why there is negative water pressure, what is the reason,
2. how to avoid the occurrence of negative water pressure.
3. Is negative water pressure due to a problem with my script.

And There is my code:
##################
from yade import ymport, utils, plot
import math
import random
from pylab import *

print '\nRunning!\n'
PACK='10Kspheres' 
intR=1.245       
DFN='penny_R0.1'
DENS=2000  
YOUNG=20e9 
ALPHA=0.12 
TENS=40e6  
COH=40e6  
FRICT=30   
KFluid=2.2e9      
visc=1.e-3                  
pFactor=1e-22   
slotAperture=1e-5 
DENS_FLUID=1000   

flowRate=0.001 

saveData=10  
iterMax=10   
saveVTK=10   
OUT=PACK+'_PlaneTests' 

O.bodies.append(ymport.text('try.txt',color=(.9,.8,.6)))

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
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
Rmean=R/numSpheres

print 'Pre-processing concluded!\n'
print 'X=',X,' | Y=',Y,' | Z=',Z,' \nSphere numbers =',numSpheres,' | Rmean=',Rmean

O.reset() 

def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,tensileStrength=TENS,cohesion=COH,frictionAngle=radians(FRICT),
                 jointNormalStiffness=YOUNG/(pi*Rmean),jointShearStiffness=ALPHA*YOUNG/(pi*Rmean),jointTensileStrength=0.,jointCohesion=0.,jointFrictionAngle=radians(FRICT),jointDilationAngle=radians(0))
def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))

mn,mx=Vector3(xinf+0.1*Rmean,yinf+0.1*Rmean,zinf+0.1*Rmean),Vector3(xsup-0.1*Rmean,ysup-0.1*Rmean,zsup-0.1*Rmean)

walls=utils.aabbWalls(oversizeFactor=1.5,extrema=(mn,mx),thickness=0.1*min(X,Y,Z),material=wallMat)
wallIds=O.bodies.append(walls)

O.bodies.append(ymport.text('try.txt',material=sphereMat,color=(.9,.8,.6)))

import gts
v1 = gts.Vertex(0,0.1,0.8)
v2 = gts.Vertex(0,0.9,0.8)
v3 = gts.Vertex(0.9,0.9,0.2)
v4 = gts.Vertex(0.9,0.1,0.2)

e1 = gts.Edge(v1,v2)
e2 = gts.Edge(v2,v4)
e3 = gts.Edge(v4,v1)
f1 = gts.Face(e1,e2,e3)

e4 = gts.Edge(v4,v3)
e5 = gts.Edge(v3,v2)
f2 = gts.Face(e2,e4,e5)

s1 = gts.Surface()
s1.add(f1)
s1.add(f2)

facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)

yade.qt.View()
yade.qt.Controller()

execfile('identifBis.py')

flow=DFNFlowEngine(    
        isActivated=False
        ,useSolver=3  
        ,boundaryUseMaxMin=[0,0,0,0,0,0] 
        ,bndCondIsPressure = [0,0,0,0,0,0]
        ,bndCondValue=[0,0,0,0,0,0] 
        ,permeabilityFactor=pFactor
        ,viscosity=visc
        ,fluidBulkModulus=KFluid
        ,clampKValues=False
        ,jointsResidualAperture=slotAperture        
)
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,label='interactionPhys')],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,neverErase=1,recordCracks=True,Key='recordcrack',label='interactionLaw')]
    ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep()),
        flow,
    PyRunner(iterPeriod=50,initRun=True,command='saveFlowVTK()',label='saveFlow',dead=1),
        NewtonIntegrator(damping=0.4,label="newton"),
        
]

def saveFlowVTK():
 flow.saveVtk(folder='injection')
 print('Save flow ', O.iter)



O.step()

SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

flow.isActivated=1
saveFlow.dead=0
O.step()
flow.imposeFlux((0.5,0.5,0.5),-flowRate)
O.dt=1e-7
########################

Best regards and thanks in advance


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.