New question #630910 on Yade: https://answers.launchpad.net/yade/+question/630910
Hi everybody, I am trying to determine which cohesive contact law I should be using in order to model rock fragmentation. So I have been implementing tests on JcfPM, CohFrictMat and CpmMat. I have adapted the uniaxial tension-Compression Test from CpmMat examples in order to use it with both CohfrictMat and JcfPM. Seems to work not too badly with JcfPM (although I can't see fracture), but there is an issue with the CohfrictMat. When I try to run the test there is an initial value of averageStress that is very high, up to -102754.2 pa. This values quickly goes back to around zero, and then increases in an 'expectable' way, but that messes up the curve a lot, and I haven't been able to find to what it is due ? I have tried setting the equilibrium distance following this thread : https://answers.launchpad.net/yade/+question/266828, but that didn't change much ( plus I wouldn't expect initial overlapping forces to be that great ?) I have also noticed that another script posted by a Yade user, also an adaptation of the CPM uniaxialcompression test but using CohFrictMat, on this thread https://answers.launchpad.net/yade/+question/372295, may also be facing the same problem. The averageStress that I get running this code ( modifying the packing in order for the warning about porosity to disappear) is -1027113730.45 Pa. Then this value also decreases inexpectantly. The adaptation with JCFpm shows no sign of such a problem, although I am using the exact same set of parameters.... Any ideas on what is going on ? Here is the code # -*- 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}') 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, }) # default parameters or from table readParamsFromTable(noTableOk=True, # unknownOk=True, young=1e9, poisson=0.25, densityCohFrictMat = 2600., FrictAngSphere = 30.*pi/180., NCohesionCohFrictMat = 4000, SCohesionCohFrictMat = 4000, intRadius=1.5, dtSafety=.8, damping=0.4, strainRateTension=.0005, strainRateCompression=.005, setSpeeds=True, # 1=tension, 2=compression (ANDed; 3=both) doModes=3, specimenLength=0.05, specimenRadius=0.1, sphereRadius=3.5e-3, # isotropic confinement (should be negative) isoPrestress=0. ) from yade.params.table import * if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description'] # make geom; the dimensions are hard-coded here; could be in param table if desired O.materials.append( CohFrictMat( young=young, frictionAngle=FrictAngSphere, poisson=poisson, density=densityCohFrictMat, normalCohesion=NCohesionCohFrictMat, shearCohesion=SCohesionCohFrictMat, fragile = True, label = 'concreteId')) sps=SpherePack() sp=pack.randomDensePack( pack.inCylinder( (0,0,-.5*specimenLength),(0,0,.5*specimenLength),specimenRadius), spheresInCell=2000,radius=sphereRadius, memoizeDb='/tmp/ triaxPackCache.sqlite',returnSpherePack=True) sp.toSimulation(material='concreteId') bb=uniaxialTestFeatures(axis = 2) 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) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb')],verletDist=.05*sphereRadius), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=intRadius,label='ss2sc')], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True)], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm = True, label='cohlaw')] ), NewtonIntegrator(damping=damping,label='damper'), UniaxialStrainer ( strainRate=strainRateTension,axis=axis,asymmetry=0, posIds=posIds,negIds=negIds, crossSectionArea=crossSectionArea,blockDisplacements=False, blockRotations=False,setSpeeds=setSpeeds,label='strainer'), PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True), PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'), ] for b in O.bodies : b.dynamic=False O.step() for i in O.interactions : i.phys.unp = i.geom.penetrationDepth for b in O.bodies : b.dynamic=True yade.qt.Controller(), yade.qt.View() #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 plot.plots={'eps':('sigma',)} #,'sigma.50')},'t':('eps')} #'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 plot.plot(subPlots=False) #O.run() initTest() waitIfBatch() -- 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 : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
