uniax.py 7.27 KB
Newer Older
1
# -*- encoding=utf-8 -*-
2
from __future__ import division
3
from __future__ import print_function
4

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

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



16
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
"""
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
44
readParamsFromTable(noTableOk=True, # unknownOk=True,
45
46
47
48
49
50
	young=24e9,
	poisson=.2,

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

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

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

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

69
70
from yade.params.table import *

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


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

79
80
81
82
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)
83
bb=uniaxialTestFeatures()
84
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
85
O.dt=dtSafety*PWaveTimeStep()
86
print('Timestep',O.dt)
87

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

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

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

111
O.saveTmp('initial');
112

113
O.timingEnabled=False
114

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

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

138
139
140
141
	O.run()

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