Uniaxial test with CohFrictMat
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:/
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:/
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.
#matplotlib.
def addPlotData():
yade.plot.addData ({'t':O.
'sigma.
'sigma.
'sigma.
})
# default parameters or from table
readParamsFromT
young=1e9,
poisson=0.25,
densityCohFrictMat = 2600.,
FrictAngSphere = 30.*pi/180.,
NCohesionCohFr
SCohesionCohFr
intRadius=1.5,
dtSafety=.8,
damping=0.4,
strainRateTens
strainRateComp
setSpeeds=True,
# 1=tension, 2=compression (ANDed; 3=both)
doModes=3,
specimenLength
specimenRadius
sphereRadius=
# isotropic confinement (should be negative)
isoPrestress=0.
)
from yade.params.table import *
if 'description' in O.tags.keys(): O.tags[
# make geom; the dimensions are hard-coded here; could be in param table if desired
O.materials.append(
CohFrictMat(
young=young,
frictionAngle=
poisson=poisson,
density=
normalCohesion
shearCohesion=
fragile = True,
label = 'concreteId'))
sps=SpherePack()
sp=pack.
pack.inCylinder(
(0,0,-
spheresInCell=
memoizeDb='/tmp/ triaxPackCache.
sp.toSimulation
bb=uniaxialTest
negIds,
O.dt=dtSafety*
print 'Timestep',O.dt
mm,mx=[pt[axis] for pt in aabbExtrema()]
coord_25,
area_25,
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
NewtonIntegrat
UniaxialStrainer
(
PyRunner(
PyRunner(
]
for b in O.bodies :
b.dynamic=False
O.step()
for i in O.interactions :
i.phys.unp = i.geom.
for b in O.bodies :
b.dynamic=True
yade.qt.
#O.miscParams=
# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that includes confinement, though
plot.plots=
O.saveTmp(
O.timingEnabled
global mode
mode='tension' if doModes & 1 else 'compression'
def initTest():
global mode
print "init"
if O.iter>0:
O.wait();
O.loadTmp(
print "Reversing plot data"; plot.reverseData()
else: plot.plot(
strainer.
try:
from yade import qt
renderer=
renderer.
except ImportError: pass
print "init done, will now run."
O.step(); # to create initial contacts
# now reset the interaction radius and go ahead
ss2sc.
is2aabb.
#O.run()
def stopIfDamaged():
global mode
if O.iter<2 or not plot.data.
sigma,
extremum=
minMaxRatio=0.5 if mode=='tension' else 0.5
if extremum==0: return
import sys; sys.stdout.flush()
if abs(sigma[
if mode=='tension' and doModes & 2: # only if compression is enabled
mode=
O.save(
print "Saved /tmp/uniax-
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.
return
else:
print "Damaged, stopping."
ft,fc=
print 'Strengths fc=%g, ft=%g, |fc/ft|
title=
print 'gnuplot'
print 'Bye.'
O.pause()
#sys.exit(0) # results in some threading exception
plot.plot(
#O.run()
initTest()
waitIfBatch()
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Jan Stránský
- Solved:
- Last query:
- Last reply: