Wrong permeability values during compaction

Hi,
The example script oedometer.py is doing exactly what you are asking and apparently you re-used it.
But at the same time you seem to have doubts on walls numbering, why? What did you change in the script?
B

Hi Bruno,

Thank you so much for your kind reply.

You are right. I am using the oedometer.py

I used its packing part and the #B part which uses FolwEngine for obtaining the permeability.

Then, I removed other parts and just used this simple vertical load: top.state.vel = (0, 0, -0.06) for the compaction purpose.

As you can see in the packing part, the walls are defined as:

walls=aabbWalls([mn,mx],thickness=0,material='walls')

wallIds=O.bodies.append(walls)

and for defining the vertical load, I used:

# walls
walls = aabbWalls()
top = walls[5]
O.bodies.

So, the top wall (vertical load) is number 5, and Qin and Qout are 2 and 3.

Would you please let me know which boundary condition and wall number I should choose in my modeling?

Best regards,
Soheil

I don't know what you "should" use. It all depends on what you want to simulate. Which direction for the pressure gradient? Which stress in each direction?

Independently: I see, now, that you are inserting walls twice, that's the main change compared to original script. I really don't understand why doing that. Prescribing boundary conditions just needs to assign triax.goal's and stressMask as you wish, not to add more walls.

B

Dear Bruno,

Thank you very much for your kind reply.

In the oedometer.py after initial triaxial compaction, in the oedometer section there is also triaxial load and unload for compaction.

But I need just vertical load for compaction. Because of this I add the ‘top.state.vel = (0, 0, -0.06)’

To be honest I am not sure which direction I should select for the pressure gradient to obtain proper results for the permeability!!!

I do not know if I should select the top and down walls or side walls!!!

I would really appreciate it if you could help me regarding this issue.

Would you please let me know what the unit of triax.goal is?

And what is stressMask and its unit?

Best regards,
Soheil

> Would you please let me know what the unit of triax.goal is?

See [1] and links therein.
TLTR: Yade itself is unit-less. Use basic SI units, it is a safe appcoach.

> And what is stressMask and its unit?

stressMask [2] is unit-less, just 0 or 1.
For convenience and compaction, packed to one integer

Cheers
Jan

[1] https://answers.launchpad.net/yade/+question/699997
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TriaxialStressController.stressMask

> in the oedometer section there is also triaxial load and unload for compaction.
> But I need just vertical load for compaction. Because of this I add the ‘top.state.vel = (0, 0, -0.06)’

An "oedometer" implies a uniaxial deformation, that's what the script does.
And you want a uniaxial deformation.
I don't know why you are trying something else.

B

Dear Jan,

Thank you very much for your kind reply and sharing the links.

They were really helpful.

Cheers,
Soheil

Dear Bruno,

Thank you very much for your kind reply.

I fixed the wall number problem and now it works properly.

But I still get weird values for the permeability.

I do not know what the issue is. Because the permeability should have a really small number like the oedometer.py example. I could gain the small value just in the first iteration and after that it increased surprisingly.

# iter permeability porosity time
1000 0.01918070551830009 0.40189088549692115 2.554930150811923
3000 564.9408214756033 0.3889401788504288 2.754930150812345
5000 371.23771416004183 0.3754329114171155 2.954930150812767
7000 418.8671964612367 0.3616597792495562 3.1549301508131893
9000 502.25383359539313 0.34705109690924985 3.3549301508136113
11000 613.7106945556807 0.33216990357515797 3.5549301508140334
13000 679.398330939237 0.31659321147545794 3.7549301508144555
15000 582.7989699519044 0.29873633244629394 3.9549301508148775
17000 601.5434724622095 0.2823106604884174 4.154930150814612
19000 627.2095474981151 0.26721453576039933 4.354930150814146

Here is my code:

###

from __future__ import print_function
from builtins import range
from yade import pack, plot
from yade import export

num_spheres=400 # number of spheres
young=1e6
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(1,1,2) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=1e99,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0.1,material='walls')
wallIds=O.bodies.append(walls)
top = walls[5]
bottom = walls[4]
leftb = walls[3]
rightb = walls[1]
rightf = walls[2]

psdSizes,psdCumm=[.01,0.03,0.04,.05,.06,.08,.12],[0.,0.1,0.3,0.3,.3,.7,1.]
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

yade.qt.views()

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.2)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1, label="flow"), #introduced as a dead engine for the moment, see 2nd section
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001:
    break

setContactFriction(radians(finalFricDegree))

triax.dead = True # (!)
top.state.vel = (0, 0, -0.3)
bottom.state.vel = leftb.state.vel = rightb.state.vel = rightf.state.vel =(0,0,0) # zero velocity of all other walls

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead = 0
flow.defTolerance = 0.3
flow.meshUpdateInterval = 200
flow.useSolver = 3
flow.permeabilityFactor = 1
flow.viscosity = 10
flow.bndCondIsPressure = [0, 0, 1, 1, 0, 0]
flow.bndCondValue = [0, 0, 1, 0, 0, 0]
flow.boundaryUseMaxMin = [0, 0, 0, 0, 0, 0]
O.dt = 0.1e-3
O.dynDt = False

O.engines=O.engines+[PyRunner(iterPeriod=2000,command='flow.saveVtk()')]

# export spheres function every 5000 iterations
O.engines += [PyRunner(command="exportSpheres()", iterPeriod=2000, initRun=True)]

def exportSpheres():
    i = O.iter
    fName = f"yade-spheres-{i:06d}.txt"
    export.text(fName)

def get_permeability():
    Qin = flow.getBoundaryFlux(5)
    Qout = flow.getBoundaryFlux(4)
    permeability = abs(Qin) / 1.e-4 #size is one, we compute K=V/∇H
    print("Qin=", Qin, " Qout=", Qout, " permeability=", permeability)
    return permeability

# porosity
O.engines += [PyRunner(iterPeriod=2000,command="plotAddData()",initRun=True)] # runs every 5000 iterations

def plotAddData():
    porosity = utils.porosity()
    permeability = get_permeability()
    plot.addData(iter=O.iter, time=O.time, porosity=porosity, permeability=permeability)
    plot.saveDataTxt("result.txt")

# run
O.stopAtIter = 20100
O.run()

###

Many thank in advance

Best regards,
Soheil

If I understand well, you are compressing a packing with an imposed velocity of 0.3.
As a result, the fluid goes out through the boundaries where pressure is imposed (top and bottom).
You then take one of these outward flux and you deduce permeability from it. But in fact the flux in your case is controlled by the imposed velocity (0.3). That's not how permeability can be determined.

For correct evaluation see the example script you started from and alternate steps B and C.
You need to make sure no deformation is going on when you evaluate permeability. A sure way is to fix DOFs and set all velocities to zero.

Bruno

Dear Bruno,

Thank you very much for your kind reply and valuable information.

Unfortunately, I could not catch this part of your answer:“You need to make sure no deformation is going on when you evaluate permeability. A sure way is to fix DOFs and set all velocities to zero.”

How can I do this? In which part?

I tried to follow the oedometer.py and changed the boundary flux, so I set the:

Qin = flow.getBoundaryFlux(2) , Qout = flow.getBoundaryFlux(3)

and I add the #C part:

#C. now the oedometer test, drained at the top, impermeable at the bottom plate
flow.bndCondIsPressure = [0, 0, 0, 1, 0, 0]
flow.bndCondValue = [0, 0, 0, 0, 0, 0]
flow.updateTriangulation = True #force remeshing to reflect new BC immediately
newton.damping = 0

Now I have received the permeability in e-11 order which is reasonable. But still I have received alternating values for the permeability. In some steps it increased and in others it decreased but the porosity decreased always.

Best regards,
Soheil

Dear Bruno,

I tried to fix the DOFs and set all velocities to zero in the permeability definition but still I give alternative values !!!

###
def get_permeability():
    for b in O.bodies:
        b.state.blockedDOFs = "xyzXYZ"
        b.state.vel = (0,0,0)
        Qin = flow.getBoundaryFlux(2)
        Qout = flow.getBoundaryFlux(3)
        permeability = abs(Qin) / 1.e-4 #size is one, we compute K=V/∇H
        print("Qin=", Qin, " Qout=", Qout, " permeability=", permeability)
        return permeability
###

I am grateful for any hints or insights into this issue.

Many thanks in advance.

Best regards,
Soheil

In the above code fragment (#11) you block DOFs but you don't solve any fluid flow after that, so it's completely ineffective.

In #10 you set "flow.bndCondIsPressure = [0, 0, 0, 1, 0, 0]", which means there is an inlet but no outlet. I would expect 0-flux from this.

You should step back and check detailed results. For instance, do you get Qin = Qout? else why? Etc.
So iteratively you will understand what goes wrong.

Bruno

Dear Bruno,

Thank you very much for your kind reply and useful suggestions.

I will go through them.

Best regards,
Soheil