yade-users team mailing list archive

[Roshan <question695717@xxxxxxxxxxxxxxxxxxxxx>]

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.