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 and set dt=0 when create the interaction as # 630910. 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

#######code####
isoPrestress=0
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.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(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=20000
]
#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[axis]-dim[axis]
for b in O.bodies:
if isinstance(b.shape,Sphere):
if (b.state.pos[axis]<(dim[axis]+layerSize*height)) or (b.state.pos[axis]> (dim[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):

def stopIfDamaged():
if O.iter<2 or 'sigma' not in plot.data : return
sigma=plot.data['sigma']
extremum=max(sigma)
if abs(strainer.strain)>5e2 : ## stop the test if the sigma is too small or strain is too large ##abs(sigma[-1]/extremum)<0.5 or
print("damaged")
O.pause()
#####count crack number
crackCount=0
for i in O.interactions:
if i.phys.cohesionBroken==True:
crackCount+=1
print('crack number=',crackCount)
plot.plots={'eps':('sigma'),'eps ':('crackNumber')}
plot.plot()
O.saveTmp()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2020-10-19
Last query:
2020-10-19
2020-10-16 Jérôme Duriez (jduriez) said on 2020-10-12: #1

Hi,

1. Sound to me like a typical "too high loading rate" issue

2. Your crackCount is the number of currently existing interactions which are now purely frictional after having been initially cohesive-frictional (i.e. after having lost their cohesive nature because of excess forces).

I do not see a reason for this number to be monotonously increasing.

If you want to consider the past (passed history of interactions breakage) and not only the present, you can not rely on "for i in O.interactions:" since this loop only considers present interactions.

In the JCFpm model for instance, there are bodies-related variables for this purpose, like JCFpmState.nbBrokenBonds

 libby (libby-eleven) said on 2020-10-16: #2

Thanks Jerome!
1. That's true. After I decrease the strain rate, the model seems to behave normally
2. thanks for pointing out the mistake. I used 'initial bond - current bond' to count the number of cracks instead and it seems working!
############
for i in O.interactions:
initialBond+=1
print('number of initial bond=',initialBond)
currentBond=0
for i in O.interactions:
if i.phys.cohesionBroken==False:
currentBond+=1
print('crack number=',initialBond-currentBond)
#################
is that an appropriate way to track the history of interactions breakage?

3. I prepare a simple script to see the effect of how fragile works, but then error occur:' ValueError: shape mismatch: objects cannot be broadcast to a single shape'. I can't figure out why, can you help me with this ?
###########
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))
O.bodies.append([sphere((0,0,0),.5,fixed=True,material=idConcrete),sphere((0,0,-1),.5,material=idConcrete)])

### define constant variables
factor=1.5
damping=.4
strainRateTension=5
bb=uniaxialTestFeatures() ## extract the information for the uniaxial test
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=0.00005*PWaveTimeStep() ## set critical time step
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),]), ## 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)],
),
UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer'),
NewtonIntegrator(gravity=(0.,0.,-1000)),
]
O.step() ## go one step to creat interactions
bo1s.aabbEnlargeFactor=1 ## reset the interaction radius
ig2ss.interactionDetectionFactor=1
O.interactions[0,1].phys.unpMax=-1

plot.plots={'i':('f'),'i ':('e'),' i':('d')}
plot.plot()
O.saveTmp()
###############################

 Jérôme Duriez (jduriez) said on 2020-10-16: #3

2- Looks like yes, at least under the assumption you're starting from initialBond=0

3- A known issue coming from an external library, see https://gitlab.com/yade-dev/trunk/-/issues/110. You can try modifying your use of plot.plot() / abandon live plotting (?)

 libby (libby-eleven) said on 2020-10-19: #4