New question #695717 on Yade:
https://answers.launchpad.net/yade/+question/695717
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)
--
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 protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
