yade-users team mailing list archive

[loiseaurare <question630910@xxxxxxxxxxxxxxxxxxxxx>]

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.