uniax.py 7.29 KB
Newer Older
Anton Gladky's avatar
Anton Gladky committed
1 2
#!/usr/bin/python
# -*- coding: utf-8 -*-
3
from __future__ import division
4
from __future__ import print_function
5

6 7
from future import standard_library
standard_library.install_aliases()
Klaus Thoeni's avatar
Klaus Thoeni committed
8
from yade import plot,pack,timing
9 10
import time, sys, os, copy

11 12 13 14 15 16
#import matplotlib
#matplotlib.rc('text',usetex=True)
#matplotlib.rc('text.latex',preamble=r'\usepackage{concrete}\usepackage{euler}')



17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
"""
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
45
readParamsFromTable(noTableOk=True, # unknownOk=True,
46 47 48 49 50 51
	young=24e9,
	poisson=.2,

	sigmaT=3.5e6,
	frictionAngle=atan(0.8),
	epsCrackOnset=1e-4,
52
	relDuctility=30,
53 54

	intRadius=1.5,
55
	dtSafety=.8,
56
	damping=0.4,
57
	strainRateTension=.05,
58
	strainRateCompression=.5,
59 60
	setSpeeds=True,
	# 1=tension, 2=compression (ANDed; 3=both)
61
	doModes=3,
62

63
	specimenLength=.15,
64
	sphereRadius=3.5e-3,
65 66 67

	# isotropic confinement (should be negative)
	isoPrestress=0,
68 69
)

70 71
from yade.params.table import *

72
if 'description' in list(O.tags.keys()): O.tags['id']=O.tags['id']+O.tags['description']
73 74 75 76 77


# 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
78
concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))
79

80 81 82 83
sps=SpherePack()
sp=pack.randomDensePack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.25*specimenLength,.17*specimenLength),spheresInCell=2000,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)
84
bb=uniaxialTestFeatures()
85
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
86
O.dt=dtSafety*PWaveTimeStep()
87
print('Timestep',O.dt)
88

89
mm,mx=[pt[axis] for pt in aabbExtrema()]
90
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
91
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)
92 93

O.engines=[
94
	ForceResetter(),
95
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
96
	InteractionLoop(
97
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
98
		[Ip2_CpmMat_CpmMat_CpmPhys()],
99
		[Law2_ScGeom_CpmPhys_Cpm()],
100
	),
101
	NewtonIntegrator(damping=damping,label='damper'),
102
	CpmStateUpdater(realPeriod=.5),
103
	UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
104 105
	PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
	PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
106
]
107
#O.miscParams=[Gl1_CpmPhys(dmgLabel=False,colorStrain=False,epsNLabel=False,epsT=False,epsTAxes=False,normal=False,contactLine=True)]
108

109
# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that includes confinement, though
110
plot.plots={'eps':('sigma',)} #,'sigma.50')},'t':('eps')} #'sigma.25','sigma.50','sigma.75')}
111

112
O.saveTmp('initial');
113

114
O.timingEnabled=False
115

116 117 118 119 120
global mode
mode='tension' if doModes & 1 else 'compression'

def initTest():
	global mode
121
	print("init")
122
	if O.iter>0:
123
		O.wait();
124
		O.loadTmp('initial')
125
		print("Reversing plot data"); plot.reverseData()
126
	else: plot.plot(subPlots=False)
127
	strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
128 129 130
	try:
		from yade import qt
		renderer=qt.Renderer()
131
		renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
132
	except ImportError: pass
133
	print("init done, will now run.")
134
	O.step(); # to create initial contacts
135
	# now reset the interaction radius and go ahead
136
	ss2sc.interactionDetectionFactor=1.
137
	is2aabb.aabbEnlargeFactor=1.
138

139 140 141 142
	O.run()

def stopIfDamaged():
	global mode
143
	if O.iter<2 or 'sigma' not in plot.data: return # do nothing at the very beginning
144
	sigma,eps=plot.data['sigma'],plot.data['eps']
145
	extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
146 147 148
	minMaxRatio=0.5 if mode=='tension' else 0.5
	if extremum==0: return
	import sys;	sys.stdout.flush()
149
	if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
150 151
		if mode=='tension' and doModes & 2: # only if compression is enabled
			mode='compression'
152
			O.save('/tmp/uniax-tension.yade.gz')
153 154
			print("Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)")
			print("Damaged, switching to compression... "); O.pause()
155 156 157
			# 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
158
			import _thread; _thread.start_new_thread(initTest,())
159 160
			return
		else:
161
			print("Damaged, stopping.")
162
			ft,fc=max(sigma),min(sigma)
163
			if doModes==3:
164
				print('Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft)))
165
			if doModes==2:
166
				print('Compressive strength fc=%g'%(abs(fc)))
167
			if doModes==1:
168
				print('Tensile strength ft=%g'%(abs(ft)))
169
			title=O.tags['description'] if 'description' in list(O.tags.keys()) else O.tags['params']
170 171
			print('gnuplot',plot.saveGnuplot(O.tags['id'],title=title))
			print('Bye.')
172 173
			O.pause()
			#sys.exit(0) # results in some threading exception
174 175
		
def addPlotData():
176
	yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
177 178 179
		'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,
180
		})
181
plot.plot(subPlots=False)
182
#O.run()
183
initTest()
184
waitIfBatch()
185