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

closed Unloading with the WeirdTriaxControl

0 votes
29 views
asked Oct 17, 2018 by MaxWiebicke (500 points)
closed Nov 20, 2018 by MaxWiebicke

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?

closed with the note: Solved
commented Oct 17, 2018 by eudoxos (49,030 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, 2018 by MaxWiebicke (500 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, 2018 by MaxWiebicke (500 points)

1 Answer

+1 vote
answered Oct 18, 2018 by eudoxos (49,030 points)
selected Nov 20, 2018 by MaxWiebicke
 
Best answer
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.
commented Nov 20, 2018 by MaxWiebicke (500 points)

Sorry, the behaviour seems fine. My mistake.

...