Author: eudoxos
Date: 2009-08-03 21:56:37 +0200 (Mon, 03 Aug 2009)
New Revision: 1914
Added:
trunk/examples/concrete/uniaxial/
trunk/examples/concrete/uniaxial/rb.py
Modified:
trunk/py/utils.py
Log:
1. Add an example of uniaxial tension-compression test for concrete (ready long
time ago...)
Added: trunk/examples/concrete/uniaxial/rb.py
===================================================================
--- trunk/examples/concrete/uniaxial/rb.py 2009-08-03 18:53:26 UTC (rev
1913)
+++ trunk/examples/concrete/uniaxial/rb.py 2009-08-03 19:56:37 UTC (rev
1914)
@@ -0,0 +1,160 @@
+# -*- encoding=utf-8 -*-
+from __future__ import division
+
+from yade import utils,plot,pack
+import time, sys, os, copy
+
+"""
+A fairly complex script performing uniaxial tension-compression test on
hyperboloid-shaped specimen.
+
+Most parameters of the model (and of the setup) can be read from table using
yade-multi.
+
+After the simulation setup, tension loading is run and stresses are
periodically saved for plotting
+as well as checked for getting below the maximum value so far. This indicates
failure (see stopIfDamaged
+function). After failure in tension, the original setup is loaded anew and the
sense of loading reversed.
+After failure in compression, strain-stress curves are saved via
plot.saveGnuplot and we exit,
+giving some useful information like peak stresses in tension/compression.
+
+Running this script for the first time can take long time, as the specimen is
prepared using triaxial
+compression. Next time, however, an attempt is made to load
previously-generated packing
+(from /tmp/triaxPackCache.sqlite) and this expensive procedure is avoided.
+
+The specimen length can be specified, its diameter is half of the length and
skirt of the hyperboloid is
+4/5 of the width.
+
+The particle size is constant and can be specified using the sphereRadius
parameter.
+
+The 3d display has displacement scaling applied, so that the fracture looks
more spectacular. The scale
+is 1000 for tension and 100 for compression.
+
+"""
+
+
+
+# default parameters or from table
+utils.readParamsFromTable(noTableOk=True, # unknownOk=True,
+ young=24e9,
+ poisson=.2,
+ G_over_E=.20,
+
+ sigmaT=3.5e6,
+ frictionAngle=atan(0.8),
+ epsCrackOnset=1e-4,
+ relDuctility=30,
+
+ intRadius=1.5,
+ dtSafety=.8,
+ damping=0.4,
+ strainRateTension=.1,
+ strainRateCompression=1,
+ setSpeeds=True,
+ # the packing has about 50% porosity, 2×2400==4800
+ density=4800,
+ # 1=tension, 2=compression (ANDed; 3=both)
+ doModes=3,
+
+ specimenLength=.2,
+ sphereRadius=3.5e-3,
+)
+
+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
+# z-oriented hyperboloid, length 20cm, diameter 10cm, skirt 8cm
+# using spheres 7mm of diameter
+spheres=pack.triaxialPack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.25*specimenLength,.2*specimenLength),radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',young=young,poisson=poisson,frictionAngle=frictionAngle,physParamsClass='CpmMat',density=density)
+O.bodies.append(spheres)
+bb=utils.uniaxialTestFeatures()
+negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
+O.dt=dtSafety*utils.PWaveTimeStep()
+print 'Timestep',O.dt
+
+mm,mx=[pt[axis] for pt in utils.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=utils.approxSectionArea(coord_25,axis),utils.approxSectionArea(coord_50,axis),utils.approxSectionArea(coord_75,axis)
+
+O.engines=[
+ BexResetter(),
+
BoundingVolumeMetaEngine([InteractingSphere2AABB(aabbEnlargeFactor=intRadius),MetaInteractingGeometry2AABB()]),
+ InsertionSortCollider(),
+ InteractionDispatchers(
+
[ef2_Sphere_Sphere_Dem3DofGeom(distanceFactor=intRadius,label='ss2d3dg')],
+
[Ip2_CpmMat_CpmMat_CpmPhys(sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E)],
+ [Law2_Dem3DofGeom_CpmPhys_Cpm()],
+ ),
+ NewtonsDampedLaw(damping=damping,label='damper'),
+ CpmPhysDamageColorizer(realPeriod=1),
+
UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
+
PeriodicPythonRunner(virtPeriod=3e-5/strainRateTension,realLim=5,command='addPlotData()',label='plotDataCollector'),
+
PeriodicPythonRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
+]
+O.miscParams=[GLDrawCpmPhys(dmgLabel=False,colorStrain=False,epsNLabel=False,epsT=False,epsTAxes=False,normal=False,contactLine=True)]
+
+plot.plots={'eps':('sigma','sigma.25','sigma.50','sigma.75')}
+plot.maxDataLen=4000
+
+O.saveTmp('initial');
+
+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()
+ strainer['strainRate']=abs(strainRateTension) if mode=='tension' else
-abs(strainRateCompression)
+ try:
+ from yade import qt
+ renderer=qt.Renderer()
+ renderer['scaleDisplacements']=True
+ renderer['displacementScale']=(1000,1000,1000) if
mode=='tension' else (100,100,100)
+ except ImportError: pass
+ print "init done, will now run."
+ 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
+ print O.tags['id'],mode,strainer['strain'],sigma[-1]
+ import sys; sys.stdout.flush()
+ if abs(sigma[-1]/extremum)<minMaxRatio:
+ if mode=='tension' and doModes & 2: # only if compression is
enabled
+ mode='compression'
+ 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)
+
+def addPlotData():
+
yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer['strain'],'sigma':strainer['avgStress'],
+
'sigma.25':utils.forcesOnCoordPlane(coord_25,axis)[axis]/area_25,
+
'sigma.50':utils.forcesOnCoordPlane(coord_50,axis)[axis]/area_50,
+
'sigma.75':utils.forcesOnCoordPlane(coord_75,axis)[axis]/area_75,
+ })
+
+
+
+initTest()
+#while True: time.sleep(1)
+
+
+
Modified: trunk/py/utils.py
===================================================================
--- trunk/py/utils.py 2009-08-03 18:53:26 UTC (rev 1913)
+++ trunk/py/utils.py 2009-08-03 19:56:37 UTC (rev 1914)
@@ -465,8 +465,8 @@
l=procStatus('VmData'); ll=l.split(); assert(ll[2]=='kB')
return int(ll[1])
-def spheresFromFileUniaxial(filename,areaSections=10,**kw):
- """Load spheres from file, but do some additional work useful for
uniaxial test:
+def uniaxialTestFeatures(filename=None,areaSections=10,**kw):
+ """Get some data about the current packing useful for uniaxial test:
1. Find the dimensions that is the longest (uniaxial loading axis)
2. Find the minimum cross-section area of the speciment by examining
several (areaSections)
@@ -475,9 +475,13 @@
3. Find the bodies that are on the negative/positive boundary, to which
the straining condition
should be applied.
+ @param filename if given, spheres will be loaded from this file (ASCII
format); if not, current simulation will be used.
+ @param areaSection number of section that will be used to estimate
cross-section
+
Returns dictionary with keys 'negIds', 'posIds', 'axis', 'area'.
"""
- ids=spheresFromFile(filename,**kw)
+ if filename: ids=spheresFromFile(filename,**kw)
+ else: ids=[b.id for b in O.bodies]
mm,mx=aabbExtrema()
dim=aabbDim(); axis=dim.index(max(dim))
import numpy
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-...@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
_______________________________________________
yade-dev mailing list
yade-dev@lists.berlios.de
https://lists.berlios.de/mailman/listinfo/yade-dev
