DFNFlowEngine

Hello

I am investigating the effect of the presence of water on the fracture of the slope, but I have encountered the following error

negative volume for an ordinary pore (temp warning, should still be safe)
*** stack smashing detected ***: /usr/bin/python3.5 terminated
Aborted (core dumped)

######################################################
from yade import ymport, utils, plot
import math
import random
from pylab import *

#### controling parameters
packing='slopegts'
smoothContact=True

output='out'
maxIter=10000

### Fluid properties ###
KFluid=2.e9 # bulk modulus do fluid (1/compressibility)
visc=1.e-3 # viscosity of the fluid
pFactor=1.8e-11 # to scale the permeability of the rock matrix: useless if lines 133-136 are not commented (impermeable matrix) -> cf. permeametre.py: 1.8e-11 gives a permeability of 1e-16 m2 for 111_10k
slotAperture=1e-3 # initial aperture of pre-existing fracture where the injection is done
DENS_FLUID=1000 # water density

flowRate=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=packing+'_PlaneTests' # nome do arquivo a ser salvo

O.bodies.append(ymport.text(packing+'.spheres'))
## preprocessing to get dimensions of the packing
dim=utils.aabbExtrema()
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

## preprocessing to get spheres dimensions
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

O.reset()

### material definition
def sphereMat(): return JCFpmMat(type=1,density=2600,young=30e10,poisson=0.25,tensileStrength=1e6,cohesion=15.5e6,frictionAngle=radians(41.3),
     jointNormalStiffness=1e10,jointShearStiffness=1e10,jointTensileStrength=0.,jointCohesion=0.,jointFrictionAngle=radians(25),jointDilationAngle=radians(0))

## Rq: density needs to be adapted as porosity of real rock is different to granular assembly due to difference in porosity (utils.sumForces(baseBodies,(0,1,0))/(Z*X) should be equal to Gamma*g*h with h=Y, g=9.82 and Gamma=2700 kg/m3

### packing ###
O.bodies.append(ymport.text(packing+'.spheres',scale=1,shift=Vector3(0,0,0),material=sphereMat))

#### Identification of the spheres on joint (some DIY here!) -> work to do on import function textExt to directly load material properties from the ascii file
inFile=open(packing+'_90_jointedPM.spheres','r')
for line in inFile:
 if '#' in line : continue
 id = int(line.split()[0])
 onJ = int(line.split()[1])
 nj = int(line.split()[2])
 j11 = float(line.split()[3])
 j12 = float(line.split()[4])
 j13 = float(line.split()[5])
 j21 = float(line.split()[6])
 j22 = float(line.split()[7])
 j23 = float(line.split()[8])
 j31 = float(line.split()[9])
 j32 = float(line.split()[10])
 j33 = float(line.split()[11])
 O.bodies[id].state.onJoint=onJ
 O.bodies[id].state.joint=nj
 O.bodies[id].state.jointNormal1=(j11,j12,j13)
 O.bodies[id].state.jointNormal2=(j21,j22,j23)
 O.bodies[id].state.jointNormal3=(j31,j32,j33)
inFile.close

#### Boundary conditions
e=4
Xmax=0
Ymax=0
baseBodies=[]

for o in O.bodies:
 if isinstance(o.shape,Sphere):
  o.shape.color=(0.9,0.8,0.6)
  ## to fix boundary particles on ground
  if o.state.pos[2]<(4) or o.state.pos[0]<(4) or o.state.pos[0]>(111):
   o.state.blockedDOFs+='xyz'
   o.shape.color=(1,1,1)

 ## to identify indicator on top
 if o.state.pos[0]>(71) and o.state.pos[0]<(72) and o.state.pos[1]>(16) and o.state.pos[1]<(17) and o.state.pos[2]>(86) and o.state.pos[2]<(88) : #single
  refPoint=o.id
  o.shape.color=(1,0,0)
  p0=o.state.pos[0]
  p2=o.state.pos[2]

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,1,0,1,1,1] # bndCondIsPressure(=vector<bool>(6, false))
                                           # bndCondIsPressure=[left, right, bottom, top, back, front]
        # ,bndCondValue=[0,0,0,0,PRESS,0]
        ,bndCondValue=[0,0,0,10,0,0] # bndCondValue(=vector<double>(6,0))
        ,permeabilityFactor=pFactor
        ,viscosity=visc
        ,fluidBulkModulus=KFluid
        ### DFN related
        ,clampKValues=False
        ,jointsResidualAperture=slotAperture
)
#### Engines definition
interactionRadius=1. # to set initial contacts to larger neighbours
O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg')],
  [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()),
        flow,
 NewtonIntegrator(damping=.7,gravity=(0.,0,-9.81)),
        PyRunner(iterPeriod=100,initRun=True,command='crackCheck()',label='check'),
        PyRunner(iterPeriod=100,initRun=True,command='saveFlowVTK()',label='saveFlow',dead=1),
        PyRunner(iterPeriod=100,initRun=True,command='saveAperture()',label='saveAperture',dead=1),
 VTKRecorder(iterPeriod=500,initRun=True,fileName=OUT+'-',recorders=['spheres','bstresses','cracks'],Key=OUT,label='saveSolid',dead=0),
 PyRunner(iterPeriod=100,initRun=False,command='jointStrengthDegradation()'),
 PyRunner(iterPeriod=10,initRun=True,command='displacement()'),
]

def crackCheck():
 flow.updateTriangulation=True

def saveFlowVTK():
 flow.saveVtk(folder='flow')

from yade import export
vtkExporter = export.VTKExporter('cracks')
def saveAperture():
 vtkExporter.exportContactPoints(what=[('b','i.phys.isBroken'),('n','i.geom.normal'),('s','i.phys.crossSection'),('a','i.phys.crackJointAperture')])

#### displacement
f = open("displacement.txt", "w")
f.write('O.iter O.time Vhorizental Vvertical Xhorizental Zvertical'+ '\n')
def displacement():

 x=O.bodies[refPoint].state.vel[0],O.bodies[refPoint].state.vel[2],O.bodies[refPoint].state.pos[0]-p0,O.bodies[refPoint].state.pos[2]-p2

 str8=str(O.iter)+' '+str(O.time)+' '+str(x)+' '
 f.write(str8+'\n')

#### joint strength degradation
stableIter=2000
stableVel=0.001
degrade=True
def jointStrengthDegradation():
 for i in O.interactions:
  if i.phys.isOnJoint :
   if i.phys.isCohesive:
    i.phys.isCohesive=False
    i.phys.FnMax=0.
    i.phys.FsMax=0.

# 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)
ss2d3dg.interactionDetectionFactor=-1.
is2aabb.aabbEnlargeFactor=-1.

### hydraulic loading
flow.isActivated=1

O.step()
#flow.imposeFlux((0.5,0.5,0.5),-flowRate)

#### RUN!!!
O.run(maxIter)