# confused by UniaxialStrainer's parameters and uniaxial compression test

Dear all,
As a beginner of Yade,I want to start with the uniaxial compression test.I refer to the code of many people on the forum.
When I run the following code, no matter how to modify the parameters of UniaxialStrainer(), I can't get the ideal stress-strain curve.The problems are as follows:
1.When the following code is run, the strain increases and the stress remains constant after an instantaneous increase.
2.When I changed the parameter asymmetry of uniaxialstrain() to 1, the cylinder was scattered from below, which puzzled me.
3.In the above process, I am also puzzled by the positive and negative values of stress and strain.In my understanding, when strainrate is set to a negative value, the stress and strain are both compressive, then both are negative. But the actual situation doesn't seem to be like this. I don't know if I understand it right.
Can anyone help me? Thanks!!

#######code####
normalCohesion=1e9
shearCohesion=1e9
RollingStiffness=1

material_1=O.materials.append(CohFrictMat(young=23e9,poisson=0.4,frictionAngle=atan(0.5),density=2036,normalCohesion=normalCohesion,shearCohesion=shearCohesion,fragile=True))
sp.toSimulation(material=material_1)

for b in O.bodies:
b.shape.color=(0,1,0)

factor=1.5
damping=0.4
strainRate_Tension=-0.005

renderer=qt.Renderer()
renderer.dispScale=(10,10,10)

bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=factor,label='ig2ss')],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False)],
),
NewtonIntegrator(damping=damping,label='damper'),
PyRunner(iterPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]

O.dt=0
O.step()
bo1s.aabbEnlargeFactor=1
ig2ss.interactionDetectionFactor=1
O.dt=0.8*PWaveTimeStep()

for i in O.interactions:
i.phys.cohesionDisablesFriction=True

O.engines=O.engines[:4]+[UniaxialStrainer(strainRate=strainRate_Tension,axis=axis,asymmetry=1,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer')]+O.engines[4:]

dim=utils.aabbExtrema()
if strainRate_Tension<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):

def stopIfDamaged():
if O.iter<2 or 'sigma' not in plot.data : return
sigma=plot.data['sigma']
extremum=max(sigma)
if extremum==0 : return
if abs(sigma[-1]/extremum)<0.5 or abs(strainer.strain)>5e2:
print('damaged')
O.pause()
plot.plots={'eps':('sigma'),'i':('eps')}
plot.plot()
O.saveTmp()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Ziyu Wang
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2021-07-28: #1

Hello,

welcome :-)

TL;TR: use "for i in O.interactions: i.phys.unp = i.geom.penetrationDepth" [1] trick (see below)

> no matter how to modify the parameters of UniaxialStrainer()

Have a look at stressUpdateInterval parameter [2] (although it is not so important).
Set stressUpdateInterval=1 to get smooth curve.
Set a larger value for larger simulations (stress evaluation needs some time, usually it not necessary to do the evaluation every iteration, but rather e.g. only before plot.addData)

> 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)

Do not use this, use posIds and negIds directly if you have them already:
for i in posIds + negIds:
O.bodies[i].shape.color = (1,1,1)

> I can't get the ideal stress-strain curve.

How do you define "ideal stress-strain curve"?

> 1.When the following code is run, the strain increases and the stress remains constant after an instantaneous increase.

The problem is that even with zero strain rate, your packing is not balanced, with high initial stress.
The next straining is too weak compared to this initial stress.
Repair it with [1]:
###
# somewhere after O.step(), before running
for i in O.interactions:
i.phys.unp = i.geom.penetrationDepth
###
The "unp trick" was discussed several times in another questions.
Doing this, the packing is in a stress-free state.

> 2.When I changed the parameter asymmetry of uniaxialstrain() to 1, the cylinder was scattered from below, which puzzled me.

By default, asymmetry is 0, i.e. the sample is loaded by symetric displacement of both ends.
For asymmetry=1, one end is fixed and the other one end is displaced.

> 3.In the above process, I am also puzzled by the positive and negative values of stress and strain.In my understanding, when strainrate is set to a negative value, the stress and strain are both compressive, then both are negative. But the actual situation doesn't seem to be like this.

"doesn't seem to be like this" - please try to make "evidence" instead of "seem".
E.g., printing strainer.strain and strainer.avgStress, both values are negative.
Maybe the plot, made by the command plot.addData(...,eps=-strainer.strain), confuses you (note the minus sign)?

Note about stopIfDamaged function, try to print the values used, the condition IMO does something else than is should.

Cheers
Jan

 Revision history for this message Jérôme Duriez (jduriez) said on 2021-08-09: #2

Regarding the sign convention of UniaxialStrainer, it seems to be as usual, conforming the continuum mechanics convention (which you exposed). See the following MWE:

#### UniaxialStrainer for compressing a 2 spheres example
for z in [0,1.99]:
O.bodies.append(sphere((0,0,z),1,dynamic=False))
O.dt = 1.e-3
O.dynDt = False
O.engines = O.engines + [UniaxialStrainer(strainRate = -1,negIds=[0],posIds= [1],crossSectionArea=1,blockDisplacements=True,label='uas')]

O.run(100,True)

signStrR = 'positive' if uas.strainRate > 0 else 'negative'
lengthVar = 'increasing' if O.bodies[uas.posIds[0]].state.pos[2]-O.bodies[uas.negIds[0]].state.pos[1] > uas.originalLength else 'decreasing'
print('With a',signStrR,'strain rate, length is',lengthVar,'and stress - strain are',uas.avgStress,uas.strain)
##########

With Yade 2021-06-21.git-c1be24c, I do get all negative values for that compressive example:

"With a negative strain rate, length is decreasing and stress - strain are -1871100.0000000698 -0.09900000000000386"

 Revision history for this message Ziyu Wang (ziyuwang1) said on 2021-08-15: #3

Your reply has helped me solve the above problems.With the solution you provided, I solved problems 1 and 2, and also figured out the problem of sign.Thank you again!
However, when I ran the modified code, I found that the stress-strain curve of uniaxial test remained linear, and there was no sign of failure.I guess it's due to the material parameters.I hope someone can help me.
Thanks!

Below is my script.

#######code#######
isoPrestress=0
normalCohesion=1e9
shearCohesion=1e9
RollingStiffness=1

material_1=O.materials.append(CohFrictMat(young=1e9,poisson=0.4,frictionAngle=atan(0.5),density=2036,normalCohesion=normalCohesion,shearCohesion=shearCohesion,fragile=False))
sp.toSimulation(material=material_1)

for b in O.bodies:
b.shape.color=(0,1,0)

factor=1.5
damping=0.4
strainRate_Tension=-0.05

#renderer=qt.Renderer()
#renderer.dispScale=(10,10,10)

bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=factor,label='ig2ss')],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False)],
),
NewtonIntegrator(damping=damping,label='damper'),
PyRunner(iterPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]

O.dt=0
O.step()

for i in O.interactions:
i.phys.unp=i.geom.penetrationDepth

bo1s.aabbEnlargeFactor=1
ig2ss.interactionDetectionFactor=1
O.dt=0.5*PWaveTimeStep()

for i in O.interactions:
i.phys.cohesionDisablesFriction=True

O.engines=O.engines[:4]+[UniaxialStrainer(strainRate=strainRate_Tension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer',stressUpdateInterval=1)]+O.engines[4:]

dim=utils.aabbExtrema()

for i in posIds + negIds:
O.bodies[i].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 extremum==0 : return
if abs(sigma[-1]/extremum)<0.3 or abs(strainer.strain)>5e2:
print('Max stress:',extremum)
print('damaged')
O.pause()

plot.setLiveForceAlwaysUpdate=True
plot.plots={'eps':('sigma'),'i':('eps')}
plot.plot()
O.saveTmp()

 Revision history for this message Jan Stránský (honzik) said on 2021-08-15: #4

> I found that the stress-strain curve of uniaxial test remained linear, and there was no sign of failure.
> I guess it's due to the material parameters.
>
> young=1e9
> normalCohesion=1e9
> shearCohesion=1e9
> RollingStiffness=1

I also guess is due to material parameters.
Have you tried to change them?

I would start with "extremes", e.g. zero cohesion, rolling stiffness etc. should give non-cohesive "weak" material which definitely should fail.
Then you can start to increase the values until you get desired response.

For experimenting, you can also increase strain rate of the loading. Maybe the specimen fail much later than you compute.

Cheers
Jan

 Revision history for this message Ziyu Wang (ziyuwang1) said on 2021-08-17: #5

Thanks Jan Stránský, that solved my question.