New question #693413 on Yade:
https://answers.launchpad.net/yade/+question/693413
Dear all,
I've been trying to model the cohesive frictional model under uniaxial force.
The problems are:
1. when the normalCohesion is low (eg.1e7), the stress-strain relationship
behaves wired, going up and down rapidly. At first, I thought it was due to the
repulsive force at the beginning, so I make i.phys.unp =
i.geom.penetrationDepth as #266828[1] and set dt=0 when create the interaction
as # 630910[2]. but still cannot resolve this issue.
2. when I set both normalCohesion and shearCohesion to 1e9 and start counting
the crack number, I find that the number of crack fluctuate a lot. this made me
confused because I've set the setCohesionOnNewContacts=False, how would the
number of cracks decrease if there was no new bond forming? By the way, the
pattern of the crack is also a bit strange, what I am expecting is that near
the peak strength, there supposed to be a surge in the crack number and remain
steady after the peak.
3. also the model seems have no difference whether I set fragile as True or
False, any suggestions is appreciated!
Thanks,
Libby
[1].https://answers.launchpad.net/yade/+question/266828
[2].https://answers.launchpad.net/yade/+question/630910
###code
from yade import pack,qt,plot
isoPrestress=0
sphereRadius=.015
normalCohesion=1e9
shearCohesion=1e9
RollingStiffness=1
idConcrete=O.materials.append(CohFrictMat(young=23.1e9,poisson=0.4,frictionAngle=atan(0.5),density=2036,normalCohesion=normalCohesion,shearCohesion=shearCohesion,fragile=False))
spp=randomDensePack(inAlignedBox((0,0,0),(.5,.5,2)),radius=sphereRadius,rRelFuzz=1/3,returnSpherePack=True,spheresInCell=2000,material=idConcrete)
# memoizeDb='/tmp/bending2DPackCache.sqlite'
spp.toSimulation()
for b in O.bodies:
b.shape.color=(0,1,0)
### define constant variables
factor=1.5
damping=.4
strainRateTension=50
## render the scale of the displacement
renderer=qt.Renderer()
renderer.dispScale=(10,10,10) ## set the render scale
bb=uniaxialTestFeatures() ## extract the information for the uniaxial test
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
#O.dt=0.5*PWaveTimeStep() ## set critical time step
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),],verletDist=.05*sphereRadius),
## expand the collision detection to creat initial interactions
##
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=factor,label='ig2ss'),
## create collision geometry
],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
## creat collision phys
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False)],
),
#GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
NewtonIntegrator(damping=damping,label='damper'),
#UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer'),##
apply strain on the both end of the body along the axis and the strain speed
is set at the beginning directly
PyRunner(iterPeriod=1,command='addPlotData()'),
PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'), ##
set the stop criteria for the simulation
#qt.SnapshotEngine(fileBase='3d',iterPeriod=1000,label='snapshot'),
#PyRunner(command='finish()',iterPeriod=2
]
#sigma=strainer.avgStress+isoPrestress
## set equlibrum for overlapping spheres
#for b in O.bodies:
#b.dynamic=False
O.dt=0
O.step() ## go one step to creat interactions
bo1s.aabbEnlargeFactor=1 ## reset the interaction radius
ig2ss.interactionDetectionFactor=1
O.dt=0.3*PWaveTimeStep()
for i in O.interactions:
i.phys.cohesionDisablesFriction=True
i.phys.unp = i.geom.penetrationDepth
O.engines=O.engines[:4]+[UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer')]+O.engines[4:]
#for b in O.bodies:
#b.dynamic=True
### reinforcing the boundary
dim=utils.aabbExtrema()
if strainRateTension>0:
layerSize=.05
height=dim[1][axis]-dim[0][axis]
for b in O.bodies:
if isinstance(b.shape,Sphere):
if (b.state.pos[axis]<(dim[0][axis]+layerSize*height)) or
(b.state.pos[axis]> (dim[1][axis]-layerSize*height)):
b.shape.color=(1,1,1)
for i in O.interactions:
if isinstance(O.bodies[i.id1].shape,Sphere) and
isinstance(O.bodies[i.id2].shape,Sphere):
if O.bodies[i.id1].shape.color==(1,1,1) or
O.bodies[i.id2].shape.color==(1,1,1):
i.phys.normalAdhesion*=1000
i.phys.shearAdhesion*=1000
def stopIfDamaged():
if O.iter<2 or 'sigma' not in plot.data : return
sigma=pl
