yade-users team mailing list archive

[Guilherme das Neves Seguro <question691351@xxxxxxxxxxxxxxxxxxxxx>]

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.