yade-users team mailing list archive

[Anay Joshi <question664499@xxxxxxxxxxxxxxxxxxxxx>]

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.