yade-users team mailing list archive

[Jan Stránský <question286479@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Jan Stránský proposed the following answer:
Hi Seti,
you have to define epsCrackOnset parameter of CpmMat (and maybe also sigmaT)
cheers
Jan


2016-02-22 12:13 GMT+01:00 Seti <question286479@xxxxxxxxxxxxxxxxxxxxx>:

> Question #286479 on Yade changed:
> https://answers.launchpad.net/yade/+question/286479
>
> Description changed to:
> Hi all,
>
> I have updated CPM model in order to simulate breaking clumps as per below
> however I face error when run the script the error is
>  FATAL /build/buildd/yade-1.07.0/pkg/dem/ConcretePM.cpp:403 go:
> Verification `!isnan(epsFracture)' failed!
> FATAL /build/buildd/yade-1.07.0/pkg/dem/ConcretePM.cpp:403 go: in
> interaction #6+#388
> terminate called without an active exception
> Aborted (core dumped)"
>
> Would you please assist?
>
> from __future__ import division
>
> from yade import plot,pack,timing
> import time, sys, os, copy
> # 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=.05,
>         strainRateCompression=.5,
>         setSpeeds=True,
>         # 1=tension, 2=compression (ANDed; 3=both)
>         doModes=2,
>
>         specimenLength=.15,
>         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']
>
> O.materials.append(CpmMat(young=30e9,poisson=0.3,frictionAngle=radians(30),density=2600,label='spheres'))
>
> O.materials.append(CpmMat(young=30e9,poisson=0.3,frictionAngle=0,density=0,label='walls'))
>
>
> walls=utils.aabbWalls([Vector3(-0.085,-0.085,-0.085),Vector3(0.085,0.085,0.085)],thickness=0.01,oversizeFactor=1.5)
> wallIds=O.bodies.append(walls)
>
>
> sp=pack.SpherePack()
> psdSizes=[0.002,0.003,0.004,0.005,0.006,0.007,0.008] # (sizes or radii of
> the grains vary from 2mm to 9.5mm)
> psdCumm=[0.14,0.28,0.34,0.50,0.65,0.85,1.00] # for the code not using
> percentage, e.g. yade-daily
>  sp.makeCloud(Vector3(-0.085,-0.085,-0.085),Vector3(0.085,0.085,0.085),-1,0,1000,False,
> 0.95,psdSizes,psdCumm,False,seed=1)
> sp.toSimulation(material='spheres')
> 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)
>
> O.engines=[
>         ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
>         InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
>                 [Ip2_CpmMat_CpmMat_CpmPhys()],
>                 [Law2_ScGeom_CpmPhys_Cpm()],
>         ),
>         NewtonIntegrator(damping=damping,label='damper'),
>         CpmStateUpdater(realPeriod=.5),
>
> 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'),
> ]
> 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)
>
> O.engines=[
>         ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
>         InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
>                 [Ip2_CpmMat_CpmMat_CpmPhys()],
>                 [Law2_ScGeom_CpmPhys_Cpm()],
>         ),
>         NewtonIntegrator(damping=damping,label='damper'),
>         CpmStateUpdater(realPeriod=.5),
>
> 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'),
> ]
>
> #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
>
> 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)
> #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     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp
>

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.