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,
