New question #702407 on Yade:
https://answers.launchpad.net/yade/+question/702407
Hi,
I am trying to do HM coupling with cylinder sample.Under the confining pressure
of 30MPa,I have get the result of which the flow pressure is 10,20,30MPa.
However,when I add the flow pressure to 40 and 50MPa,I get the following error
prompt and my simulation stop:
###
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
ThreadRunner:35 void yade::ThreadRunner::run(): Exception
occured:
CGAL ERROR: assertion violation!
Expr: dexp != 2047
File: /usr/include/CGAL/Mpzf.h
Line: 411
Explanation: Creating an Mpzf from infinity or NaN.
###
I do not know what happened.I will attach my code in the end.Thanks for anyone
help!
###code###
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys
damp=0.4
young=66.2e9
name='JCFPM_triax'
poisson=0.522
name='cylinder'
preStress=-30e6
strainRate = -0.1
OUT=str(preStress)+'_JCFPM_triax'
r1=0.002
r2=0.01
nw=24
nh=15
rParticle=0.000731723
bcCoeff = 5
width = 0.025
height = 0.05
allx,ally,allz,r=np.genfromtxt('packing-cylinder.spheres', unpack=True)
mnx=min(allx)*1.01 #same as aabbExtrema, get the xyz of the lowest corner,
multiply by factor to reduce it to let the wall surround all the spheres.
mny=min(ally)*1.01
mnz=min(allz)*0.99
mxx=max(allx)*1.01 #same as aabbExtrema, get the xyz of the top corner,
multiply by factor to increase it to let the wall surround all the spheres.
mxy=max(ally)*1.01
mxz=max(allz)*1.01
Sy=50e6
O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=8e6,cohesion=40e6,frictionAngle=radians(25),label='sphere'))
O.materials.append(JCFpmMat(type=1,label='wall1',young=young,poisson=poisson,frictionAngle=radians(25)))
O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='wall2'))
mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0,material='wall2')
wallIds=O.bodies.append(walls)
bottom,top=Vector3(0,0,0),Vector3(0,0,0.05)
sp=pack.SpherePack()
pred=pack.inCylinder(bottom,top,radius)
sp=pack.randomDensePack(pred,radius=0.0005,rRelFuzz=0,returnSpherePack=True,memoizeDb='/tmp/triax.sqlite')
spheres=sp.toSimulation(color=(0.9,0.8,0.6),material='sphere')
bot = [O.bodies[s] for s in spheres if
O.bodies[s].state.pos[2]height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
for s in bot:
s.shape.color = (1,0,0)
s.state.blockedDOFs = 'xyzXYZ'
#s.state.vel = (0,0,-vel)
for s in top:
s.shape.color = Vector3(0,1,0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0,0,vel)
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
for h in range(nh):
v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),
rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),
rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),
rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),
rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
f1 = facet((v1,v2,v3),color=(0,0,1),material='wall1')
f2 = facet((v1,v3,v4),color=(0,0,1),material='wall1')
facets.extend((f1,f2))
O.bodies.append(facets)
mass = O.bodies[7].state.mass
for f in facets:
f.state.mass = mass
f.state.blockedDOFs = 'XYZz'
def lateral():
elatTot=0.0
nTot=0
for b in O.bodies:
x=b.state.refPos[0]
y=b.state.refPos[1]
d=sqrt(pow(x,2)+pow(y,2))
if d > r1 and d < r2:
b.shape.color=(0.6,0.5,0.15)
xnew=b.state.pos[0]
ynew=b.state.pos[1]
dnew=sqrt(pow(xnew,2)+pow(ynew,2))
elat=(dnew-d)/d
elatTot+=elat
nTot+=1
elat_avg=elatTot/nTot
return elat_avg
def addForces():
for f in facets:
n = f.shape.normal
a = f.shape.area
O.forces.addF(f.id,preStress*a*n)
def stopIfDamaged(maxEps=0.001):
extremum = max(abs(s) for s in plot.data['s'])
s = abs(plot.data['s'][-1])
e = abs(plot.data['e'][-1])
if O.iter < 1000 or e < maxEps:
return
if abs(s)/abs(extremum) < 0.5 :
print('Simulation finished')
presentcohesive_count = 0
for i in O.interactions:
if hasattr(i.phys, 'isCohesive'):
if i.phys.isCohesive == True:
presentcohesive_count+=1
print('the number of cohesive bond now
is:',presentcohesive_count)
print('Max stress and strain:',extremum,e)
O.pause()
O.dt=.5*utils.PWaveTimeStep()
newton=NewtonIntegrator(d
