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:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-10-16
Last reply:

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.