yade-users team mailing list archive

[Chareyre <question691351@xxxxxxxxxxxxxxxxxxxxx>]

Question #691351 on Yade changed:
https://answers.launchpad.net/yade/+question/691351

    Status: Open => Answered

Chareyre proposed the following answer:
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 <
question691351@xxxxxxxxxxxxxxxxxxxxx> 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     : yade-users@xxxxxxxxxxxxxxxxxxx
> 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>

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