# Why doesn't the strain equal to zero

Asked by Zhicheng Gao on 2021-06-07

I use the O.cell.trsf=Matrix3.Identity to set the current cell configuration to be the reference one, but when I use print(O.cell.trsf, triax.strain) to print the current strain, the strain matrix is not equal to zero, this is my code:

this is my code:
##______________ First section, generate sample_________

from __future__ import print_function
from yade import pack, qt, plot
from math import *

nRead=readParamsFromTable(
## model parameters
num_spheres=100,
targetPorosity= .387,
confiningPressure=-100000,
## material parameters
compFricDegree=15,#contact friction during the confining phase
finalFricDegree=30,#contact friction during the deviatoric loading
young=2e8,
poisson=.2,
density=2600,
alphaKr=7.5,
alphaKtw=0,
competaRoll=.22,
finaletaRoll=.22,
etaTwist=0,
normalCohesion=0,
shearCohesion=0,
## fluid parameters
fluidDensity=1000,
dynamicViscosity=.001,
## control parameters
damp=0,
stabilityThreshold=.001,
## output specifications
filename='suffusion',
unknowOk=True
)

from yade.params.table import *

O.periodic=True
O.cell.hSize=Matrix3(.001,0,0, 0,.001,0, 0,0,.001)
# create materials for spheres
#shear strength is the sum of friction and adhesion, so the momentRotationLaw=True
O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres'))

# generate particles packing
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.001,.001,.001),psdSizes=[0.00008,0.000125,0.0001592,0.0002003,0.0003153,0.000399,0.000502,0.0005743],psdCumm=[0.0,0.00628,0.0565,0.198,0.721,0.915,0.991,1.0],num=num_spheres,seed=1)
sp.toSimulation(material='spheres')

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)],
),
PeriodicFlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
PeriTriaxController(label='triax',
# specify target values and whether they are strains or stresses
goal=(confiningPressure,confiningPressure,confiningPressure), stressMask=7,
# type of servo-control, the strain rate isn't determined, it shloud check the unbalanced force
dynCell=True,maxStrainRate=(10,10,10),
# wait until the unbalanced force goes below this value
maxUnbalanced=stabilityThreshold,relStressTol=1e-3,
doneHook='compactionFinished()'
),
NewtonIntegrator(damping=0)
]
qt.View()

import sys
def compactionFinished():
#check the current porosity
# if the current porosity is lager than target Porosity and comFricDegree is lager than 10,
# then we decrease friction value and apply it to all the bodies and contacts,
# else we decrease rolling friction value.
global compFricDegree, competaRoll
if porosity()>targetPorosity and compFricDegree>5:
# we decrease friction value and apply it to all the bodies and contacts
compFricDegree = 0.95*compFricDegree
setContactFriction(radians(compFricDegree))
print('Friction:', compFricDegree,'porosity:', porosity())
# python syntax, make each step printout
sys.stdout.flush()
elif porosity()>targetPorosity:
# we decrease rolling fiction value and apply it to all the bodies and contacts
competaRoll=0.95*competaRoll
for b in O.bodies:
b.mat.etaRoll=competaRoll
for i in O.interactions:
i.phys.etaRoll=competaRoll
print('Rollingfriction:', b.mat.etaRoll, 'porosity:', porosity())
sys.stdout.flush()
else:
# after sample preparation, save the state
O.save('compactedState'+filename+'.yade.gz')
print('Compacted state saved', 'porosity', porosity())
# next time, called python command
triax.doneHook=''
O.pause()

# enable energy tracking in the code
O.trackEnergy=True

# define function to record history
def history():
plot.addData(unbalanced=unbalancedForce(),i=O.iter,exx=-triax.strain,
eyy=-triax.strain, ezz=-triax.strain,
sxx=-triax.stress,syy=-triax.stress,szz=-triax.stress,
ev=-triax.strain-triax.strain-triax.strain,
porosity=porosity(),Etot=O.energy.total(),**O.energy# save all available energy data
)
O.engines=O.engines+[PyRunner(command='history()', iterPeriod=20)]
# define what to plot
plot.plots={'i':('unbalanced','porosity'),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),' i ':('Etot')}
# show the plot
plot.plot()
plot.live=True

O.run()
O.wait()

plot.saveDataTxt('compactedState'+filename)
##__________________________________second section, deviatoric loading__________________________
# change the contact parameters to the final calibration value
setContactFriction(radians(finalFricDegree))
for b in O.bodies:
b.mat.etaRoll=finaletaRoll
for i in O.interactions:
i.phys.etaRoll=finaletaRoll
print(O.cell.hSize,O.cell.trsf)
# set the current cell configuration to be the reference one
O.cell.trsf=Matrix3.Identity
print(O.cell.trsf, triax.strain)
# change control type: keep constant confinement in x,y, 30% compression in z
triax.goal=(confiningPressure,confiningPressure,-.3)
triax.stressMask=3
# allow faster deformation along x,y to better maintain stresses
triax.maxStrainRate=(1.,1.,.01)
# call triaxFinished instead of compactionFinished
triax.doneHook='triaxFinished()'

def triaxFinished():
O.save('loadedState'+filename+'.yade.gz')
print('Finished')
O.pause()

# Reset all plot data; keep plots and labels intact
plot.resetData()
O.run(1000,True)
plot.saveDataTxt('loadedState'+filename)

## 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:
 Revision history for this message Jan Stránský (honzik) said on 2021-06-07: #1

Hello,
please read  and try to provide code in the form of MWE

M = minimal.
This code is far from minimal, vast majority of the code is not related to the problem at all - most of parameters, materials, PyRunner functions, FlowEngine...
A few particles, triax, a few steps, is enough. Then illustrating your problem on specific command.

W = working.
I copy-pasted your code, started it and got:
if porosity()>targetPorosity and compFricDegree>5:
IndentationError: unexpected indent

Cheers
Jan

 Revision history for this message Zhicheng Gao (zhichenggao) said on 2021-06-07: #2

I'm sorry to provide a lot of useless code, the key code is ：
O.cell.trsf=Matrix3.Identity
print(O.cell.trsf, triax.strain)
Then, I found that the deformation matrix is the identity matrix, but the triaxial strain is not zero.

 Revision history for this message  Jan Stránský (honzik) said on 2021-06-07: #3

a) "modifying the position directly is likely to break Yade’s algorithms" . O.cell.trsf is similar as setting position directly.
b) Some values are "propagated" to correct values during the iteration.

PeriTriaxController.strain is only updated inside the engine's action [3,4]. So you can:
a) try to run one more step to let triax be updated
b) set cell.velGrad to get trsf=I in the next step

Cheers
Jan

 Revision history for this message Zhicheng Gao (zhichenggao) said on 2021-06-07: #4

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

To post a message you must log in.