Why is there so much negative water pressure

Asked by Liu Changdong on 2020-10-16

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

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-10-27
Last reply:
2020-10-28

Hi,
Did you check that a flux equal to -0.01 doesn't actually means fluid is pumped out (instead of being injected)?
Bruno

Liu Changdong (changdong) said : #2

Hi Bruno
Thank you very much for your answer. But I don't quite understand your answer. Isn't [1] injecting flux at postion (x,y,z)? And positive for outflux , negative for influx . Can you explain that again?

tkanks
Changdong Liu

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.imposeFlux

I honestly don't remember the sign convention but I guess it is easy to check on your side.
If changing the sign of the source term changes the sign of the fluid pressure then you have you answer.
If it doesn't, then something else is going on.
If you find inconsistent sign convention between documentation and code please report it here.
Bruno

Can you help with this problem?

Provide an answer of your own, or ask Liu Changdong for more information if necessary.

To post a message you must log in.