Permeability error

Hi,
I can't make sense of problem outline.

> there is a plane of fractures identified and the flow only occurs in this plane direction (y axis)

One plane with multiple fractures? How do you define orientation of a plane by one single vector (y), do you mean the unit normal?

> the flow will occurs due to a pressure difference imposed on the lower facet (y_min)

Mmmh... if y is the normal of the fracture should gradient be along y too?

> acceptable interval of difference between Qin and Qout

It is unclear what you call Qin and Qout.

> I've put all files together on [1]

We need one single file, not four, and no external link, please [1].
Regards
Bruno

[1] https://www.yade-dem.org/wiki/Howtoask

Hello mr Chareyre;

> I can't make sense of problem outline.

It's a cubic rock matrix modelled by DEM with a fracture plane in it. There will be a pressure difference imposed on its bottom so I can measure its permeability.

> One plane with multiple fractures? How do you define orientation of a plane by one single vector (y), do you mean the unit normal?

Yes. The plane is defined by reading the file "surfacePlane" I've mentioned. The surface is generated by a Matlab code.

> Mmmh... if y is the normal of the fracture should gradient be along y too?

Yes.

> It is unclear what you call Qin and Qout.
Qin and Qout are the flow rates calculated at the inflow and outflow boundaries respectively.

> We need one single file, not four, and no external link, please [1].
I've put all files together because the main code calls the other ones and it's also too long. I'll put the main code down here anyways.

Thanks for all the help!

---------------------------------------

# -*- coding: utf-8 -*-
# encoding: utf-8
from yade import ymport, utils, plot
import math
import random
from pylab import *

print '\nRunning!\n'

### Packages e DFN (Discrete Fracture Network)
PACK='10Kspheres' # spheres
intR=1.245 # to define the average number of bonds per particle -> cf. coordination number.py
DFN='penny_R0.1'

### parameters of the simulation (geometry, boundary conditions, output, etc...)

### Material microproperties (set of parameters for simulating Colton sandstone with K=10)
DENS=4000
YOUNG=30e9
ALPHA=0.33
TENS=40e5
COH=40e6
FRICT=40
PRESS=0.0

### Mechanical loading (state of stress before the injection) ###
Sxx=-4.9e6 # Sigmaxx
Syy=-3.9e6 # Sigmayy
Szz=-4.1e6 # Sigmazz

### Fluid properties ###
KFluid=2.e9 # bulk modulus do fluidd (1/compressibility)
visc=1.e-3 # viscosity of the fluid
pFactor=1.8e-11
slotAperture=1e-3 # initial aperture of pre-existing fracture where the injection is done
DENS_FLUID=1000 # water density

### hydraulic loading ###
flowRate=0 # injection flow rate (original 1e-5)
bottom=0

### Simulation Control ###
saveData=10 # data record interval
iterMax=10 # numero maximo de iteracoes da simulacao (passos de tempo)
saveVTK=10 # number of Vtk files
OUT=PACK+'_PlaneTests' # nome do arquivo a ser salvo
# print 'Carregado: pacote de esferas, tensoes, propriedades e flowrate!\n'

### Pre-processing ###
O.bodies.append(ymport.text(PACK+'.spheres'))

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

###
# all previous lines were for getting dimensions of the packing to create walls at the right positions (below) because walls have to be generated after spheres for FlowEngine
O.reset()
###

#### here we reconstruct the scene with right dimensions (because walls have to be imported before spheres for flow engine)

### material definition
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))

### walls ###
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)

### packing ###
O.bodies.append(ymport.text(PACK+'.spheres',material=sphereMat))

################ DFN ################
print '\nReading surface...\n'
arq = open("surfacePlane.txt", "r")

# Read squares
linha = arq.readline()
valores = linha.split()
squares = int(valores[0])

# set coordinates
xxx = np.zeros((4,1),dtype = np.float)
yyy = np.zeros((4,1),dtype = np.float)
zzz = np.zeros((4,1),dtype = np.float)
center=np.zeros((squares,3),dtype = np.float)

# create surface
surf = gts.Surface()
i=0
total=0
for linha in arq:
    valores = linha.split()
    xxx[i]=(float(valores[0]))
    yyy[i]=(float(valores[1]))
    zzz[i]=(float(valores[2]))
    i=i+1
    if (mod(i,4) == 0):
        center[total,0]=(xxx[0]+xxx[1]+xxx[2]+xxx[3])/4
        center[total,1]=(yyy[0]+yyy[1]+yyy[2]+yyy[3])/4
        center[total,2]=(zzz[0]+zzz[1]+zzz[2]+zzz[3])/4
        total = total + 1
        i = 0
        pt1 = gts.Vertex(xxx[0], yyy[0], zzz[0])
        pt2 = gts.Vertex(xxx[1], yyy[1], zzz[1])
        pt3 = gts.Vertex(xxx[2], yyy[2], zzz[2])
        pt4 = gts.Vertex(xxx[3], yyy[3], zzz[3])
        e1 = gts.Edge(pt1,pt2)
        e2 = gts.Edge(pt2,pt3)
        e3 = gts.Edge(pt3,pt1)
        face = gts.Face(e1,e2,e3)
        surf.add(face)
        e1 = gts.Edge(pt1,pt3)
        e2 = gts.Edge(pt3,pt4)
        e3 = gts.Edge(pt4,pt1)
        face = gts.Face(e1,e2,e3)
        surf.add(face)
print '\nTotal faces: ',total
arq.close()
surface=gtsSurface2Facets(surf,wire=False,material=wallMat)
O.bodies.append(surface)
print '\nCalls "identifyInitialFractures.py" \n(identification of interaction on pre-existing fractures)\n'
execfile('identifyInitialFractures.py') # call to identification of interaction on pre-existing fractures

### ENGINES ###

### Triaxial Engine ###
triax=TriaxialStressController(
 internalCompaction=False
 ,stressMask=7
 ,goal1=Sxx
 ,goal2=Syy
 ,goal3=Szz
 ,max_vel=0.01
)

### Flow Engine ###
flow=DFNFlowEngine(
        isActivated=False
        ,useSolver=3 # (0: Gauss Seidel, 1: Taucs, 2: Pardiso, 3: CHOLMOD)
        #,boundaryUseMaxMin=[0,0,0,0,0,0] # [left, right, bottom, top, back, front]: False means boundary made with walls
        ,bndCondIsPressure = [0,0,bottom,bottom,0,0] # bndCondIsPressure(=vector<bool>(6, false))
                                           # bndCondIsPressure=[left, right, bottom, top, back, front]
        # ,bndCondValue=[0,0,0,0,PRESS,0]
        ,bndCondValue=[0,0,PRESS,0,0,0] # bndCondValue(=vector<double>(6,0))
        ,permeabilityFactor=pFactor
        ,viscosity=visc
        ,fluidBulkModulus=KFluid
        ### DFN related
        ,clampKValues=False
        ,jointsResidualAperture=slotAperture
)

# With DFNFlow, we can block every cells not concerned with fractures with the following function:
# Rk: if these lines are commented (permeable matrix), we get warning about cholmod: is it an issue?
# I am not sure yet but it is annoying and it seems to slow the calculation...
def blockStuff():
 for k in range(flow.nCells()): flow.blockCell(k,True)
flow.blockHook="blockStuff()"

### Definicao da simulacao ###
# (DEM loop, interaction law, servo control, recording, etc...)
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=OUT,label='interactionLaw')]
 ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep()),
        triax,
        flow,
        NewtonIntegrator(damping=0.4,label="newton"),
        PyRunner(iterPeriod=int(1),initRun=True,command='crackCheck()',label='check'),
        # PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='recData'),
        PyRunner(iterPeriod=int(1),initRun=True,command='saveFlowVTK()',label='saveFlow',dead=1),
        PyRunner(iterPeriod=int(1),initRun=True,command='saveAperture()',label='saveAperture',dead=1),
        VTKRecorder(iterPeriod=int(1),initRun=True,fileName=OUT+'-',recorders=['spheres','bstresses','cracks'],Key=OUT,label='saveSolid',dead=0)
]

### custom functions ###
# check if new cracks are created to update "flow mesh permeability":
# Rk: if commented (with associated pyrunner), results are totally different???
cks=cks0=0
def crackCheck():
  global tensCks, shearCks, cks, cks0
  cks=interactionLaw.nbTensCracks+interactionLaw.nbShearCracks
  if cks>(cks0):
    # print 'new crack! Update triangulation!'
    flow.updateTriangulation=True
  cks0=cks

### save flow field (pressure and velocity)
def saveFlowVTK():
 flow.saveVtk(folder='injection')
 #print('Salvou flow ', O.iter)

### save cracks aperture
from yade import export
vtkExporter = export.VTKExporter('cracks')
def saveAperture():
  print 'Saved aperture in ', O.iter
  vtkExporter.exportContactPoints(what=[('b','i.phys.isBroken'),('n','i.geom.normal'),('s','i.phys.crossSection'),('a','i.phys.crackJointAperture')])

### save macroscopic data
ex0=ey0=ez0=0
def recorder():
    global ex0,ey0,ez0
    crackVolume=crackSurface=0
    for i in O.interactions:
        if i.phys.isBroken:
     crackVolume+=i.phys.crossSection*i.phys.crackJointAperture
     crackSurface+=i.phys.crossSection
    yade.plot.addData( t=O.time
   ,i=O.iter
      # Deformacoes
   ,ex=triax.strain[0]-ex0
   ,ey=triax.strain[1]-ey0
   ,ez=triax.strain[2]-ez0
      # Tensoes principais aplicadas no bloco
   ,sx=triax.stress(triax.wall_right_id)[0]
   ,sy=triax.stress(triax.wall_top_id)[1]
   ,sz=triax.stress(triax.wall_front_id)[2]
   ,p=flow.getPorePressure((xinf+X/2.,yinf+Y/2.,zinf+Z/2.)) # Mexer aqui - colocar as coordenadas centrais de um quadrado do plano de fratura
   ,tc=interactionLaw.nbTensCracks
   ,sc=interactionLaw.nbShearCracks
   ,p32=crackSurface
   ,p33=crackVolume
   ,unbF=utils.unbalancedForce() # unbalancedforce a ser carregada mais a frente
    )
    plot.saveDataTxt(OUT)
    print ('\nSaved macroscopic data in ', O.iter)

########## SIMULATION #########

# Simulation starts here
### manage interaction detection factor during the first timestep (near neighbour bonds are created at first timestep)
O.step()
## initializes the interaction detection factor to default value (new contacts, frictional, between strictly contacting particles)
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.
saveSolid.dead=1

### mechanical loading
while 1:
  O.run(100, True)
  print 'unbalanced force = ',unbalancedForce()
  if ( unbalancedForce()<0.005 and ((abs(abs(triax.stress(triax.wall_right_id)[0])-abs(Sxx))/abs(Sxx))<0.001) and ((abs(abs(triax.stress(triax.wall_top_id)[1])-abs(Syy))/abs(Syy))<0.001) and ((abs(abs(triax.stress(triax.wall_front_id)[2])-abs(Szz))/abs(Szz))<0.001) ):
    print '\nStabilizing || iteration=', O.iter
    O.run(100,True) # to further stabilize the system
    print '\nConfined state \nSxx=',triax.stress(triax.wall_right_id)[0],' | Syy=',triax.stress(triax.wall_top_id)[1],' | Szz=',triax.stress(triax.wall_front_id)[2]
    ex0=triax.strain[0]
    ey0=triax.strain[1]
    ez0=triax.strain[2]
    O.save(OUT+'_confined.yade')
    break

### hydraulic loading
print 'Activate flow engine now || iteration = ', O.iter
triax.max_vel=1
flow.isActivated=1
saveFlow.dead=0
saveSolid.dead=0
saveAperture.dead=1
O.step() # needed to avoid segfault?
saveFlow.iterPeriod=int(iterMax/saveVTK)
saveSolid.iterPeriod=int(iterMax/saveVTK)
saveAperture.iterPeriod=int(iterMax/saveVTK)

### Fluid injection on pre-existing fractures ###
print '\nApply fluid on fractures'
# imposeFlux((Vector3)pos, (float)p) → None
cont =0
for i in range(0, squares):
    flow.imposeFlux((center[i,0],center[i,1],center[i,2]),-flowRate)
    # Impose a flux in cell located at ‘pos’ (i.e. add a source term in the flow problem).
    # Outflux positive, influx negative.
    cont=cont+1

print '\nFlow rate = ', flowRate, '\nInjection in ', cont, ' squares\n'

plot.plots={'i':('p',None,'tc')}

print '\n1st iteractions\n'
iter = 1
while iter <=iterMax:
 O.run(int(1),True)
 # print O.iter
 iter=iter+1

O.run(int(iterMax), True)

print '\n2nd iteractions'
print '\nRedefine Flow Engine\n'
bottom = 1
PRESS = 3e6

O.run(1,True)
# getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(2)
Qout = flow.getBoundaryFlux(3)
permeability = abs(Qout)*flow.viscosity*Y/(X*Z) # !!! if Pout=1, Pin=0
conductivity = permeability*DENS_FLUID*9.82/flow.viscosity # K=rho*g*k/nu
print "\nQin=",Qin,"\nQout=",Qout,"\nOBS: ARE THEY EQUAL? IF NOT => NO FLOW!\n"
print "\nPermeability [m2]=",permeability,"\nHydraulic conductivity [m/s]=",conductivity, '\n\nTHE END!\n'

So the fracture plane is orthogonal to "y" and the main flow direction is
"y"?
What do you expect from this?
B

On Thu, 18 Jun 2020 at 00:51, Guilherme das Neves Seguro <
<email address hidden>> wrote:

> Question #691351 on Yade changed:
> https://answers.launchpad.net/yade/+question/691351
>
> Status: Answered => Open
>
> Guilherme das Neves Seguro is still having a problem:
> Hello mr Chareyre;
>
> > I can't make sense of problem outline.
>
> It's a cubic rock matrix modelled by DEM with a fracture plane in it.
> There will be a pressure difference imposed on its bottom so I can
> measure its permeability.
>
>
> > One plane with multiple fractures? How do you define orientation of a
> plane by one single vector (y), do you mean the unit normal?
>
> Yes. The plane is defined by reading the file "surfacePlane" I've
> mentioned. The surface is generated by a Matlab code.
>
> > Mmmh... if y is the normal of the fracture should gradient be along y
> too?
>
> Yes.
>
> > It is unclear what you call Qin and Qout.
> Qin and Qout are the flow rates calculated at the inflow and outflow
> boundaries respectively.
>
>
> > We need one single file, not four, and no external link, please [1].
> I've put all files together because the main code calls the other ones and
> it's also too long. I'll put the main code down here anyways.
>
> Thanks for all the help!
>
> ---------------------------------------
>
> # -*- coding: utf-8 -*-
> # encoding: utf-8
> from yade import ymport, utils, plot
> import math
> import random
> from pylab import *
>
> print '\nRunning!\n'
>
> ### Packages e DFN (Discrete Fracture Network)
> PACK='10Kspheres' # spheres
> intR=1.245 # to define the average number of bonds per particle ->
> cf. coordination number.py
> DFN='penny_R0.1'
>
> ### parameters of the simulation (geometry, boundary conditions, output,
> etc...)
>
> ### Material microproperties (set of parameters for simulating Colton
> sandstone with K=10)
> DENS=4000
> YOUNG=30e9
> ALPHA=0.33
> TENS=40e5
> COH=40e6
> FRICT=40
> PRESS=0.0
>
> ### Mechanical loading (state of stress before the injection) ###
> Sxx=-4.9e6 # Sigmaxx
> Syy=-3.9e6 # Sigmayy
> Szz=-4.1e6 # Sigmazz
>
> ### Fluid properties ###
> KFluid=2.e9 # bulk modulus do fluidd (1/compressibility)
> visc=1.e-3 # viscosity of the fluid
> pFactor=1.8e-11
> slotAperture=1e-3 # initial aperture of pre-existing fracture where the
> injection is done
> DENS_FLUID=1000 # water density
>
> ### hydraulic loading ###
> flowRate=0 # injection flow rate (original 1e-5)
> bottom=0
>
> ### Simulation Control ###
> saveData=10 # data record interval
> iterMax=10 # numero maximo de iteracoes da simulacao (passos de tempo)
> saveVTK=10 # number of Vtk files
> OUT=PACK+'_PlaneTests' # nome do arquivo a ser salvo
> # print 'Carregado: pacote de esferas, tensoes, propriedades e flowrate!\n'
>
> ### Pre-processing ###
> O.bodies.append(ymport.text(PACK+'.spheres'))
>
> 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
>
> ###
> # all previous lines were for getting dimensions of the packing to create
> walls at the right positions (below) because walls have to be generated
> after spheres for FlowEngine
> O.reset()
> ###
>
> #### here we reconstruct the scene with right dimensions (because walls
> have to be imported before spheres for flow engine)
>
> ### material definition
> 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))
>
> ### walls ###
>
> 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)
>
> ### packing ###
> O.bodies.append(ymport.text(PACK+'.spheres',material=sphereMat))
>
> ################ DFN ################
> print '\nReading surface...\n'
> arq = open("surfacePlane.txt", "r")
>
> # Read squares
> linha = arq.readline()
> valores = linha.split()
> squares = int(valores[0])
>
> # set coordinates
> xxx = np.zeros((4,1),dtype = np.float)
> yyy = np.zeros((4,1),dtype = np.float)
> zzz = np.zeros((4,1),dtype = np.float)
> center=np.zeros((squares,3),dtype = np.float)
>
> # create surface
> surf = gts.Surface()
> i=0
> total=0
> for linha in arq:
> valores = linha.split()
> xxx[i]=(float(valores[0]))
> yyy[i]=(float(valores[1]))
> zzz[i]=(float(valores[2]))
> i=i+1
> if (mod(i,4) == 0):
> center[total,0]=(xxx[0]+xxx[1]+xxx[2]+xxx[3])/4
> center[total,1]=(yyy[0]+yyy[1]+yyy[2]+yyy[3])/4
> center[total,2]=(zzz[0]+zzz[1]+zzz[2]+zzz[3])/4
> total = total + 1
> i = 0
> pt1 = gts.Vertex(xxx[0], yyy[0], zzz[0])
> pt2 = gts.Vertex(xxx[1], yyy[1], zzz[1])
> pt3 = gts.Vertex(xxx[2], yyy[2], zzz[2])
> pt4 = gts.Vertex(xxx[3], yyy[3], zzz[3])
> e1 = gts.Edge(pt1,pt2)
> e2 = gts.Edge(pt2,pt3)
> e3 = gts.Edge(pt3,pt1)
> face = gts.Face(e1,e2,e3)
> surf.add(face)
> e1 = gts.Edge(pt1,pt3)
> e2 = gts.Edge(pt3,pt4)
> e3 = gts.Edge(pt4,pt1)
> face = gts.Face(e1,e2,e3)
> surf.add(face)
> print '\nTotal faces: ',total
> arq.close()
> surface=gtsSurface2Facets(surf,wire=False,material=wallMat)
> O.bodies.append(surface)
> print '\nCalls "identifyInitialFractures.py" \n(identification of
> interaction on pre-existing fractures)\n'
> execfile('identifyInitialFractures.py') # call to identification of
> interaction on pre-existing fractures
>
> ### ENGINES ###
>
> ### Triaxial Engine ###
> triax=TriaxialStressController(
> internalCompaction=False
> ,stressMask=7
> ,goal1=Sxx
> ,goal2=Syy
> ,goal3=Szz
> ,max_vel=0.01
> )
>
> ### Flow Engine ###
> flow=DFNFlowEngine(
> isActivated=False
> ,useSolver=3 # (0: Gauss Seidel, 1: Taucs, 2: Pardiso, 3: CHOLMOD)
> #,boundaryUseMaxMin=[0,0,0,0,0,0] # [left, right, bottom, top,
> back, front]: False means boundary made with walls
> ,bndCondIsPressure = [0,0,bottom,bottom,0,0] #
> bndCondIsPressure(=vector<bool>(6, false))
> # bndCondIsPressure=[left,
> right, bottom, top, back, front]
> # ,bndCondValue=[0,0,0,0,PRESS,0]
> ,bndCondValue=[0,0,PRESS,0,0,0] #
> bndCondValue(=vector<double>(6,0))
> ,permeabilityFactor=pFactor
> ,viscosity=visc
> ,fluidBulkModulus=KFluid
> ### DFN related
> ,clampKValues=False
> ,jointsResidualAperture=slotAperture
> )
>
> # With DFNFlow, we can block every cells not concerned with fractures with
> the following function:
> # Rk: if these lines are commented (permeable matrix), we get warning
> about cholmod: is it an issue?
> # I am not sure yet but it is annoying and it seems to slow the
> calculation...
> def blockStuff():
> for k in range(flow.nCells()): flow.blockCell(k,True)
> flow.blockHook="blockStuff()"
>
> ### Definicao da simulacao ###
> # (DEM loop, interaction law, servo control, recording, etc...)
> 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=OUT,label='interactionLaw')]
> ),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=0.1*utils.PWaveTimeStep()),
> triax,
> flow,
> NewtonIntegrator(damping=0.4,label="newton"),
>
> PyRunner(iterPeriod=int(1),initRun=True,command='crackCheck()',label='check'),
> #
> PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='recData'),
>
> PyRunner(iterPeriod=int(1),initRun=True,command='saveFlowVTK()',label='saveFlow',dead=1),
>
> PyRunner(iterPeriod=int(1),initRun=True,command='saveAperture()',label='saveAperture',dead=1),
>
> VTKRecorder(iterPeriod=int(1),initRun=True,fileName=OUT+'-',recorders=['spheres','bstresses','cracks'],Key=OUT,label='saveSolid',dead=0)
> ]
>
> ### custom functions ###
> # check if new cracks are created to update "flow mesh permeability":
> # Rk: if commented (with associated pyrunner), results are totally
> different???
> cks=cks0=0
> def crackCheck():
> global tensCks, shearCks, cks, cks0
> cks=interactionLaw.nbTensCracks+interactionLaw.nbShearCracks
> if cks>(cks0):
> # print 'new crack! Update triangulation!'
> flow.updateTriangulation=True
> cks0=cks
>
> ### save flow field (pressure and velocity)
> def saveFlowVTK():
> flow.saveVtk(folder='injection')
> #print('Salvou flow ', O.iter)
>
> ### save cracks aperture
> from yade import export
> vtkExporter = export.VTKExporter('cracks')
> def saveAperture():
> print 'Saved aperture in ', O.iter
>
> vtkExporter.exportContactPoints(what=[('b','i.phys.isBroken'),('n','i.geom.normal'),('s','i.phys.crossSection'),('a','i.phys.crackJointAperture')])
>
> ### save macroscopic data
> ex0=ey0=ez0=0
> def recorder():
> global ex0,ey0,ez0
> crackVolume=crackSurface=0
> for i in O.interactions:
> if i.phys.isBroken:
> crackVolume+=i.phys.crossSection*i.phys.crackJointAperture
> crackSurface+=i.phys.crossSection
> yade.plot.addData( t=O.time
> ,i=O.iter
> # Deformacoes
> ,ex=triax.strain[0]-ex0
> ,ey=triax.strain[1]-ey0
> ,ez=triax.strain[2]-ez0
> # Tensoes principais aplicadas no bloco
> ,sx=triax.stress(triax.wall_right_id)[0]
> ,sy=triax.stress(triax.wall_top_id)[1]
> ,sz=triax.stress(triax.wall_front_id)[2]
>
> ,p=flow.getPorePressure((xinf+X/2.,yinf+Y/2.,zinf+Z/2.)) # Mexer aqui -
> colocar as coordenadas centrais de um quadrado do plano de fratura
> ,tc=interactionLaw.nbTensCracks
> ,sc=interactionLaw.nbShearCracks
> ,p32=crackSurface
> ,p33=crackVolume
> ,unbF=utils.unbalancedForce() # unbalancedforce a
> ser carregada mais a frente
> )
> plot.saveDataTxt(OUT)
> print ('\nSaved macroscopic data in ', O.iter)
>
> ########## SIMULATION #########
>
> # Simulation starts here
> ### manage interaction detection factor during the first timestep (near
> neighbour bonds are created at first timestep)
> O.step()
> ## initializes the interaction detection factor to default value (new
> contacts, frictional, between strictly contacting particles)
> SSgeom.interactionDetectionFactor=-1.
> Saabb.aabbEnlargeFactor=-1.
> saveSolid.dead=1
>
> ### mechanical loading
> while 1:
> O.run(100, True)
> print 'unbalanced force = ',unbalancedForce()
> if ( unbalancedForce()<0.005 and
> ((abs(abs(triax.stress(triax.wall_right_id)[0])-abs(Sxx))/abs(Sxx))<0.001)
> and
> ((abs(abs(triax.stress(triax.wall_top_id)[1])-abs(Syy))/abs(Syy))<0.001)
> and
> ((abs(abs(triax.stress(triax.wall_front_id)[2])-abs(Szz))/abs(Szz))<0.001)
> ):
> print '\nStabilizing || iteration=', O.iter
> O.run(100,True) # to further stabilize the system
> print '\nConfined state \nSxx=',triax.stress(triax.wall_right_id)[0],'
> | Syy=',triax.stress(triax.wall_top_id)[1],' |
> Szz=',triax.stress(triax.wall_front_id)[2]
> ex0=triax.strain[0]
> ey0=triax.strain[1]
> ez0=triax.strain[2]
> O.save(OUT+'_confined.yade')
> break
>
> ### hydraulic loading
> print 'Activate flow engine now || iteration = ', O.iter
> triax.max_vel=1
> flow.isActivated=1
> saveFlow.dead=0
> saveSolid.dead=0
> saveAperture.dead=1
> O.step() # needed to avoid segfault?
> saveFlow.iterPeriod=int(iterMax/saveVTK)
> saveSolid.iterPeriod=int(iterMax/saveVTK)
> saveAperture.iterPeriod=int(iterMax/saveVTK)
>
> ### Fluid injection on pre-existing fractures ###
> print '\nApply fluid on fractures'
> # imposeFlux((Vector3)pos, (float)p) → None
> cont =0
> for i in range(0, squares):
> flow.imposeFlux((center[i,0],center[i,1],center[i,2]),-flowRate)
> # Impose a flux in cell located at ‘pos’ (i.e. add a source term in
> the flow problem).
> # Outflux positive, influx negative.
> cont=cont+1
>
> print '\nFlow rate = ', flowRate, '\nInjection in ', cont, ' squares\n'
>
> plot.plots={'i':('p',None,'tc')}
>
> print '\n1st iteractions\n'
> iter = 1
> while iter <=iterMax:
> O.run(int(1),True)
> # print O.iter
> iter=iter+1
>
> O.run(int(iterMax), True)
>
> print '\n2nd iteractions'
> print '\nRedefine Flow Engine\n'
> bottom = 1
> PRESS = 3e6
>
> O.run(1,True)
> # getBoundaryFlux get the total discharge [m3/s]
> Qin = flow.getBoundaryFlux(2)
> Qout = flow.getBoundaryFlux(3)
> permeability = abs(Qout)*flow.viscosity*Y/(X*Z) # !!! if Pout=1, Pin=0
> conductivity = permeability*DENS_FLUID*9.82/flow.viscosity # K=rho*g*k/nu
> print "\nQin=",Qin,"\nQout=",Qout,"\nOBS: ARE THEY EQUAL? IF NOT => NO
> FLOW!\n"
> print "\nPermeability [m2]=",permeability,"\nHydraulic conductivity
> [m/s]=",conductivity, '\n\nTHE END!\n'
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

--
--
_______________
Bruno Chareyre
Associate Professor
ENSE³ - Grenoble INP
Lab. 3SR
BP 53
38041 Grenoble cedex 9
Tél : +33 4 56 52 86 21
________________

Email too brief?
Here's why: email charter
<https://marcuselliott.co.uk/wp-content/uploads/2017/04/emailCharter.jpg>

I think I understand the problem.

You want to measure the permeability of a fracture contained inside a "rock" sample and, for that, you want to relate the flux inside the fracture to the pressure difference along its extension. Am I right?

Let me know if I got it right:

1 - you set up a cubic sample cut by a persistent fracture plane along the Y direction
2 - you use the blockStuff() function to allow flow only in the cells that are cut by the fracture, all the other cells are thus "impermeable"
3 - you apply the pressure difference on the ymin and ymax boundaries of the cubic sample to generate flow inside the fracture only (you use the no flux condition on all the other boundaries xmin, xmax, zmin, zmax).
4 - you can see the flow along the fracture plane (confirm with the visualization of the vtk file in Paraview)
5 - when you measure Qin and Qout with the getBoundaryFlux() function on the ymin and ymax boundaries they are not equal

If all points above are right, your problem is the difference you get between Qin and Qout. I think I (i.e., someone I worked with) experienced such problem before and It may comes from the getBoundaryFlux() function and the way it works. In your case, you should have flow only in the non blocked cells contained inside the sample but I don't know on which cells the getBoundaryFlux() function measure the flux (external/virtual cells or cells inside the sample?). Bruno may help with that, otherwise you'll have to read the source code.

To see if the problem is coming from the getBoundaryFlux() function, you could check the flow velocity (or the flux) at different points located inside the fracture. If the velocity is the same everywhere, you could even use it to compute the permeability of your fracture without using the Qin and Qout.

BTW, I am not sure you can use this

abs(Qout)*flow.viscosity*Y/(X*Z)

to measure the permeability of the fracture since the flow does not go through the entire surface X*Z of the sample.

Instead, you should look into the cubic law [1] and use something like that

Q=(deltaP*Width*h**3)/(12*mu*Length)

to measure the hydraulic aperture h of the fracture before computing its permeability by

permeability=(h**2)/12

Luc

---

[1] https://inis.iaea.org/collection/NCLCollectionStore/_Public/26/034/26034259.pdf

Hello mr Scholtès;

> You want to measure the permeability of a fracture contained inside a "rock" sample and, for that, you want to relate the flux inside the fracture to the pressure difference along its extension. Am I right?

Yes, exactly. Also all 5 considerations you did were all correct.

> If all points above are right, your problem is the difference you get between Qin and Qout.

Yes. That's the trouble I'm facing and got stuck.

> I think I (i.e., someone I worked with) experienced such problem before and It may comes from the getBoundaryFlux() function and the way it works. [...] To see if the problem is coming from the getBoundaryFlux() function, you could check the flow velocity (or the flux) at different points located inside the fracture. If the velocity is the same everywhere, you could even use it to compute the permeability of your fracture without using the Qin and Qout.

About measuring the velocity on different points of the fracture should I use averageVelocity() or is there any other more suitable function for this purpose?

> Instead, you should look into the cubic law [1]

Thanks for this, I'll read the work and change it on the code. Thanks for all others advices as well, I'm gonna try all the suggestions.