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[0],
            eyy=-triax.strain[1], ezz=-triax.strain[2],
            sxx=-triax.stress[0],syy=-triax.stress[1],szz=-triax.stress[2],
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            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 :
#1

Hello,
please read [1] 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

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#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
Best Jan Stránský (honzik) said :
#3

a) "modifying the position directly is likely to break Yade’s algorithms" [2]. 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

[2] https://yade-dem.org/doc/user.html#imposed-velocity
[3] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/PeriIsoCompressor.cpp#L126
[4] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/PeriIsoCompressor.cpp#L195

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#4

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