Welcome to ask.woodem.org. You may post when you login through your GitHub account.

Unloading with the WeirdTriaxControl

0 votes
10 views
asked Oct 17 by MaxWiebicke (400 points)

Hi Vaclav,

I want to simulate load-unload cycles with the WeirdTriaxControl, but I encountered a "weird" behaviour for the unloading part. Not matter the initial conditions and parameters I set, I get a volumetric increase upon unloading, although I would expect a decrease, see the figure

I do simulate the loading - unloading with setting the stressMasks to S.lab.stressMasks=[3,7] and the goals to S.lab.goals=[-0.05,latSig]. The same happens if I do the unloading with controlling the vertical strain. In ealier simulations, three years ago when we started, I did get compressive behaviour upon unloading.

Do you have any idea what I am getting wrong?

commented Oct 17 by eudoxos (47,770 points)

Those 3 years ago, were you also using WeirdTriaxialController? Another thing(maybe a stupid question), are you really setting the engine parameters, are they set from values in S.lab.* ? If you can prepare a MWE of what you're doing, I can help perhaps better. Also, you can compile with debugging (debug=1 with scons), then run woo -D and set woo.log.setLevel('WeirdTriaxControl',woo.log.TRACE) to see copious output about what is happening.

IIRC I unified stress/strain computation from contacts across a few places (there were about 3 implementations, one of them was giving incorrect result), do you think that could explain that? In that case I would check git history to pinpoint what changed.

Sorry I can't tell you better straight away.

v.

commented Oct 18 by MaxWiebicke (400 points)

Yes, we were also using WeirdTriaxialController. And yes, I remember that you worked on unifying the stress computations -- didn't know about the strain computations, that might be an issue.
And I am using the values from S.lab to set the engine.

I produced a mwe (although it feels too long to be considered an mwe), but I didn't know how to shorten it even more.

import woo
from minieigen import *
from woo.core import *
from woo.dem import *
import sys

woo.master.usesApi=10104

S=woo.master.scene=Scene(fields=[DemField()],trackEnergy=True)
mat=FrictMat(young=50e6,ktDivKn=0.4,tanPhi=0.4)
S.plot.plots={'i':('sig_x','sig_y','sig_z'),
              'i ':('eps_x','eps_y','eps_z','epsV'), 
              'sigM':('q',),
              't':('**S.energy'),
              't ':('unbE','unbF')
              }

S.engines=DemField.minimalEngines(damping=0.2) # basic DEM engines
# load spheres from external file
sp=ShapePack()
sp.loadTxt('initial-state-reference.txt')
sp.toDem(S,S.dem,mat)
S.dem.collectNodes()

S.engines=S.engines+[PyRunner(100,'addPlotData2(S)')] # save plot data every 100 steps

S.cell.trsf=Matrix3.Identity # reset current cell strain (transformation matrix)

def addPlotData2(S):
    sig=S.lab.triax.stress.diagonal()
    unbF=S.lab.triax.currUnbalanced
    eps = S.cell.trsf.diagonal()-(1,1,1)
    S.plot.addData(sig=sig,eps=eps,i=S.step,unbF=unbF,step=S.lab.step,sigM=sig.sum()/3,q=(sig[2]-sig[0]),t=S.time,unbE=woo.utils.unbalancedEnergy(S),**S.energy,epsV=eps.sum())

epsRate=1. # deformation rate (1/s)
what=S.lab.what='triax'
# save the step we are in in the loop below
S.lab.step=-1 
S.lab._setWritable('step')

mm=sum([p.mass for p in S.dem.par]) # mass of particles, for the "mass" magic on the next line
# stressMask=3 means sig_x, sig_y (and eps_z) prescribed
S.engines=[WeirdTriaxControl(goal=(0,0,0),stressMask=3,maxStrainRate=epsRate*Vector3.Ones,absStressTol=10,relStressTol=.01,maxUnbalanced=1e-5,doneHook='nextStep(S)',label='triax',mass=100*mm)]+S.engines
S.one() # to compute initial stress, which is the prescribed lateral stress
latSig=woo.utils.stressStiffnessWork()[0].head().sum()/3.  # average stress
S.lab.triax.goal=(latSig,latSig,0)

S.lab.goals=[-0.05,latSig]
# stress mask is given by the sum of the stress values (sig_x = 1, sig_y = 2, sig_z = 4)
# 3=1+2 means control sig_x and sig_y, and 7=1+2+4 means sig_x, sig_y and sig_z 
S.lab.stressMasks=[3,7]
if (len(S.lab.stressMasks)!=len(S.lab.goals)):
    print("Different number of goals and paths given ...")
    sys.exit()

def nextStep(S):
    if S.lab.step>=len(S.lab.goals): return
    if S.lab.step>=0:
        goal=S.lab.goals[S.lab.step]
        print('%s/%d: finished eps=%g (real %g)'%(S.lab.what,S.lab.step,goal,S.cell.trsf[2,2]-1))
    S.lab.step+=1
    if S.lab.step>=len(S.lab.goals):
        S.plot.saveDataTxt('%s.plot-data.txt'%what)
        S.plot.saveGnuplot('%s.gnuplot'%what)
        S.stop()
        return
    goal=S.lab.goals[S.lab.step]
    eps_zz=S.cell.trsf[2,2]-1
    S.lab.triax.goal[2]=goal
    S.lab.triax.stressMask=S.lab.stressMasks[S.lab.step]
    print('%s/%d: true triaxial (eps %g -> %g)'%(S.lab.what,S.lab.step,eps_zz,goal))

nextStep(S)
S.run()

That poduces the following result
enter image description here

commented Oct 18 by MaxWiebicke (400 points)

1 Answer

0 votes
answered Oct 18 by eudoxos (47,770 points)
Sorry for the mistake, strain is computed still the same, but stress was updated; the relevant change is here: https://github.com/woodem/woo/commit/2ef28cafbb1610f03b327c86f98833fad2c42314 plus this is the DemFuncs::stressStiffness function which is newly used https://github.com/woodem/woo/blob/a6be2c47f918ede238ae884cba40abcac632c154/pkg/dem/Funcs.cpp#L183 -- perhaps you spot some problem in there? I will check better later, thanks for your code. v.
...