Question #664499 on Yade changed:
https://answers.launchpad.net/yade/+question/664499

Description changed to:
Hi, 
I am simulating 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

Reply via email to