DFNFlowEngine

Asked by Roshan on 2021-02-23

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)

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2021-02-23
Last reply:
2021-03-12
Launchpad Janitor (janitor) said : #1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Robert Caulk (rcaulk) said : #2

Changing status to answered in case one of the DFNFlow owners decides to answer this someday.

Can you help with this problem?

Provide an answer of your own, or ask Roshan for more information if necessary.

To post a message you must log in.