New question #664499 on Yade:
https://answers.launchpad.net/yade/+question/664499
Hi,
I am simulation a uniaxial compression of a concrete cylinder in yade. I am
Using the uniax.py reference code from the github and I am editing that as per
my problem. I have edited the code to use the JCFpm model for the cohesion. I
am copying the code and the error i am getting in the following.
CODE:
#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
from yade import plot,pack,timing
import time, sys, os, copy
#import matplotlib
#matplotlib.rc('text',usetex=True)
#matplotlib.rc('text.latex',preamble=r'\usepackage{concrete}\usepackage{euler}')
# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
young=24e9,
poisson=.2,
sigmaT=3.5e6,
frictionAngle=atan(0.8),
epsCrackOnset=1e-4,
relDuctility=30,
intRadius=1.5,
dtSafety=.8,
damping=0.4,
strainRateTension=.01,
strainRateCompression=.01,
setSpeeds=True,
# 1=tension, 2=compression (ANDed; 3=both)
doModes=2,
specimenLength=.08,
sphereRadius=1.25e-3,
# isotropic confinement (should be negative)
isoPrestress=0,
)
from yade.params.table import *
from yade import pack, plot
#############################Simulation
Control#################################
rate =-0.01 #Deformation rate
iterMax = 10000 #Maximum number of iterations
saveVTK = 2000 #Opt files for paraview
if 'description' in O.tags.keys():
O.tags['id']=O.tags['id']+O.tags['description']
rad = 5
# make geom; the dimensions are hard-coded here; could be in param table if
desired
# z-oriented hyperboloid, length 20cm, diameter 10cm, skirt 8cm
# using spheres 7mm of diameter
############################# Material Definition
##############################
concreteId=O.materials.append(JCFpmMat(type=1,young=young,frictionAngle=frictionAngle,poisson=poisson,density=2400,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress,cohesion=1e6))
sps=SpherePack()
sp=pack.randomDensePack(pack.inCylinder((0,0,1),(0,0,1.080),
0.01*rad),spheresInCell=500,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
#sp=pack.randomDensePack(pack.inAlignedBox((-.25*specimenLength,-.25*specimenLength,-.5*specimenLength),(.25*specimenLength,.25*specimenLength,.5*specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=concreteId)
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt
mm,mx=[pt[axis] for pt in aabbExtrema()]
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)
########################### Engines are defined here ###########################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.5*sphereRadius),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,
label='interactionPhys')],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True, Key=Out,
label='interactionLaw')],
),
NewtonIntegrator(damping=damping,label='damper'),
#CpmStateUpdater(realPeriod=.5),
UniaxialStrainer(strainRate=strainRateCompression,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=setSpeeds,
stopStrain=0.1, dead=1, label='strainer'),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,
defaultDt=utils.PWaveTimeStep()),
PyRunner(virtPeriod=1e-6/strainRateCompression,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]
#O.miscParams=[Gl1_CpmPhys(dmgLabel=False,colorStrain=False,epsNLabel=False,epsT=False,epsTAxes=False,normal=False,contactLine=True)]
# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that
includes confinement, though
#Recorder and Plot
def recorder():
yade.plot.addData({'i':O.iter,
'eps':strainer.strain,
'sigma':strainer.avgStress,
'tc':interactionLaw.nbTensCracks,
'sc':interactionLaw.nbShearCracks,
'te':interactionLaw.totalTensCracksE,
'se':interactionLaw.totalShearCracksE,
'unbF':utils.unbalancedForce()})
plot.saveDataTxt(OUT)
#############intraction detection#################
# interaction range ((cf. A DEM model for soft and hard rock, Scholtes &
Donze, JMPS 2013))
O.step();
#plot.plots={'eps':('sigma',)} #,'sigma.50')},'t':('eps')} Plot during the
simulation #'sigma.25','sigma.50','sigma.75')}
O.saveTmp('initial');
O.timingEnabled=False
global mode
mode='tension' if doModes & 1 else 'compression'
def initTest():
global mode
print "init"
if O.iter>0:
O.wait();
O.loadTmp('initial')
print "Reversing plot data"; plot.reverseData()
else: plot.plot(subPlots=False)
strainer.strainRate=abs(strainRateTension) if mode=='tension' else
-abs(strainRateCompression)
try:
from yade import qt
renderer=qt.Renderer()
renderer.dispScale=(1000,1000,1000) if mode=='tension' else
(100,100,100)
except ImportError: pass
print "init done, will now run."
O.step(); # to create initial contacts
# now reset the interaction radius and go ahead
ss2sc.interactionDetectionFactor=1.
is2aabb.aabbEnlargeFactor=1.
O.run()
def stopIfDamaged():
global mode
if O.iter<2 or not plot.data.has_key('sigma'): return # do nothing at
the very beginning
sigma,eps=plot.data['sigma'],plot.data['eps']
extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
minMaxRatio=0.5 if mode=='tension' else 0.5
if extremum==0: return
import sys; sys.stdout.flush()
if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if
isoPrestress==0 else 5e-2):
if mode=='tension' and doModes & 2: # only if compression is
enabled
mode='compression'
O.save('/tmp/uniax-tension.yade.gz')
print "Saved /tmp/uniax-tension.yade.gz (for use with
interaction-histogram.py and uniax-post.py)"
print "Damaged, switching to compression... "; O.pause()
# important! initTest must be launched in a separate
thread;
# otherwise O.load would wait for the iteration to
finish,
# but it would wait for initTest to return and deadlock
would result
import thread; thread.start_new_thread(initTest,())
return
else:
print "Damaged, stopping."
ft,fc=max(sigma),min(sigma)
print 'Strengths fc=%g, ft=%g,
|fc/ft|=%g'%(fc,ft,abs(fc/ft))
title=O.tags['description'] if 'description' in
O.tags.keys() else O.tags['params']
print
'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
print 'Bye.'
O.pause()
#sys.exit(0) # results in some threading exception
def addPlotData():
yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
'sigma.25':forcesOnCoordPlane(coord_25,axis)[axis]/area_25+isoPrestress,
'sigma.50':forcesOnCoordPlane(coord_50,axis)[axis]/area_50+isoPrestress,
'sigma.75':forcesOnCoordPlane(coord_75,axis)[axis]/area_75+isoPrestress,
})
plot.plot(subPlots=False)
strainer.dead=0
#O.run()
initTest()
waitIfBatch()
ERROR:
File "/apps/pkg/yade-2017.01a/rhel7_u2-x86_64/gnu/bin/yade-2017.01a", line
182, in runScript
execfile(script,globals())
File "uniax_0.5mm_jcfpm.py", line 86, in <module>
concreteId=O.materials.append(JCFpmMat(type=1,young=young,frictionAngle=frictionAngle,poisson=poisson,density=2400,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress,cohesion=1e6))
AttributeError: No such attribute: epsCrackOnset.
I am fairly a new user to this, can anyone help me understand and suggest the
solution for this?
My second question is,
Is JCFpm be the appropriate method for modeling the concrete? I planned to use
the Bonded particle model (Cundal, 2004) for the cohesion but here on this
portal, Luc Scholtes suggested that JCFpm is actually better than BPM. Can
anyone suggest me on the which method would be appropriate?
Thanks,
Anay
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
