orthogonal packing

Asked by Fedor

Hello! I try to use regular packing in the isotropic pressure during triaxial test but how to consider the young modulus in this test? I want to check how prosity will be changed during isopropic compression with orthogonal packing. But when we use cloud packing this parametr affects on result.

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
Robert Caulk (rcaulk) said :
#1

Hello,

Please try to clarify your question - it is not immediately clear what the question is. You would like to know how to record packing porosity? What does that have to do with regular packing? Why does a "cloud packing" affect "result?" What is "result?" You want to know how to measure the macroscopic young's modulus? What does the Young's modulus have to do with the porosity?

Please, and I beg you, please read this [1] if you seek assistance.

some important bullet points:
- Keep 1 question per thread
- Provide an MWE
- Copy and paste errors

Cheers,

Robert

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

Revision history for this message
Fedor (vodovozov) said :
#2

Thank you for notes. My apologies

I want to use orthogonal packing and investigate how to change porosity vs isotropic pressure with diffrent values of young modulus.
I used the existing script and added young modulus, but when I vary module the porosity changes inadequate. There is the script.
from __future__ import print_function
O.periodic = True
O.cell.hSize = Matrix3(2, 0, 0, 0, 2, 0, 0, 0, 2)

from yade import pack, plot
sigmaIso = -7e5
young=5e8
poisson=0.2
frictionAngle=0

# the "if 0:" block will be never executed, therefore the "else:" block will be
# to use cloud instead of regular packing, change to "if 1:" or something similar

 # in this case, add dense packing
pp=O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=frictionAngle))
O.bodies.append(pack.regularOrtho(pack.inAlignedBox((0, 0, 0), (2, 2, 2)), radius=.09, gap=0, color=(0, 0, 1), material=pp))
 #O.bodies.append(pack.regularHexa(pack.inAlignedBox((0, 0, 0), (2, 2, 2)), radius=.09, gap=0, color=(0, 0, 1)))

# create "dense" packing by setting friction to zero initially
#O.materials[0].frictionAngle = 0

# simulation loop (will be run at every step)
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
                # interaction loop
                [Ig2_Sphere_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(damping=.4),
        # run checkStress function (defined below) every second
        # the label is arbitrary, and is used later to refer to this engine
        PyRunner(command='checkStress()', realPeriod=1, label='checker'),
        # record data for plotting every 100 steps; addData function is defined below
        PyRunner(command='addData()', iterPeriod=100),
 PeriTriaxController(
                label='triax',
                # specify target values and whether they are strains or stresses
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=7,
                # type of servo-control
                dynCell=True,
                maxStrainRate=(10, 10, 10),
                # wait until the unbalanced force goes below this value
                maxUnbalanced=5,
                relStressTol=1e-3,
                # call this function when goal is reached and the packing is stable
                doneHook='compactionFinished()'
        ),
]

# set the integration timestep to be 1/2 of the "critical" timestep
O.dt = .5 * PWaveTimeStep()

# prescribe isotropic normal deformation (constant strain rate)
# of the periodic cell
O.cell.velGrad = Matrix3(0, 0, 0, 0, -.1, 0, 0, 0, 0)

# when to stop the isotropic compression (used inside checkStress)
limitMeanStress = -5e5

# called every second by the PyRunner engine
def checkStress():
 # stress tensor as the sum of normal and shear contributions
 # Matrix3.Zero is the intial value for sum(...)
 stress = getStress().trace() / 3.
 print('mean stress', stress)
 porrr=utils.voxelPorosity(200,(0.75,0.75,0.75),(1.25,1.25,1.25))
 #porrr=utils.aabbExtrema(0, False)
 print('porosity', porrr)

# called from the 'checker' engine periodically, during the shear phase

# called periodically to store data history
def addData():
 # get the stress tensor (as 3x3 matrix)
 stress = sum(normalShearStressTensors(), Matrix3.Zero)
 # give names to values we are interested in and save them
 plot.addData(exz=O.cell.trsf[0, 2], szz=stress[2, 2], sxz=stress[0, 2], tanPhi=(stress[0, 2] / stress[2, 2]) if stress[2, 2] != 0 else 0, i=O.iter, porrr=utils.voxelPorosity(200,(0.75,0.75,0.75),(1.25,1.25,1.25)), stress = getStress().trace() / 3.)
 # color particles based on rotation amount
 for b in O.bodies:
  # rot() gives rotation vector between reference and current position
  b.shape.color = scalarOnColorScale(b.state.rot().norm(), 0, pi / 2.)

# define what to plot (3 plots in total)
## exz(i), [left y axis, separate by None:] szz(i), sxz(i)
## szz(exz), sxz(exz)
## tanPhi(i)
# note the space in 'i ' so that it does not overwrite the 'i' entry
plot.plots = { 'i ': ('stress',), 'stress':('porrr')}

# better show rotation of particles
Gl1_Sphere.stripes = True

# open the plot on the screen
plot.plot()

O.saveTmp()

def compactionFinished():
 # set the current cell configuration to be the reference one
 O.cell.trsf = Matrix3.Identity
 # change control type: keep constant confinement in x,y, 20% compression in z
 triax.goal = (sigmaIso, sigmaIso, sigmaIso)
 triax.stressMask = 7
 # allow faster deformation along x,y to better maintain stresses
 triax.maxStrainRate = (1., 1., 1.)
 # next time, call triaxFinished instead of compactionFinished
 triax.doneHook = 'triaxFinished()'
 # do not wait for stabilization before calling triaxFinished
 triax.maxUnbalanced = 10

def triaxFinished():
 print('Finished')
 O.pause()

Revision history for this message
Jan Stránský (honzik) said :
#3

> I used the existing script

what is "the existing script"? Do you have any reference (e.g. a tutorial script for this and that page)?

> but when I vary module the porosity changes inadequate.

What does "vary" mean?
What does "inadequate" mean?
Please be (much) more specific. E.g. for this young I get this porosity, for that young I get that porosity...

Cheers
Jan

Revision history for this message
Fedor (vodovozov) said :
#4

This file is called *04-periodic-simple-shear.py*
<https://gitlab.com/yade-dev/trunk/-/blob/990dda83e2558b7d244f52476f07742ce6b44176/doc/sphinx/tutorial/04-periodic-simple-shear.py>.
The link is
https://gitlab.com/yade-dev/trunk/-/blob/990dda83e2558b7d244f52476f07742ce6b44176/doc/sphinx/tutorial/04-periodic-simple-shear.py

For module 5e7 the porosity
equals 0,39 (isotropic pressure -1e5)
equals 0,44 (isotropic pressure -5e5)
equals 0,43 (isotropic pressure -7e5)

Thank you for your help

On Tue, Nov 23, 2021 at 9:15 PM Jan Stránský <
<email address hidden>> wrote:

> Your question #699605 on Yade changed:
> https://answers.launchpad.net/yade/+question/699605
>
> Status: Open => Needs information
>
> Jan Stránský requested more information:
> > I used the existing script
>
> what is "the existing script"? Do you have any reference (e.g. a
> tutorial script for this and that page)?
>
> > but when I vary module the porosity changes inadequate.
>
> What does "vary" mean?
> What does "inadequate" mean?
> Please be (much) more specific. E.g. for this young I get this porosity,
> for that young I get that porosity...
>
> Cheers
> Jan
>
> --
> To answer this request for more information, you can either reply to
> this email or enter your reply at the following page:
> https://answers.launchpad.net/yade/+question/699605
>
> You received this question notification because you asked the question.
>

Revision history for this message
Jan Stránský (honzik) said :
#5

Thanks for clarification.

What does "inadequate" mean here in the relation of porosity vs. young?
What is expected behavior?

Cheers
Jan

Revision history for this message
Fedor (vodovozov) said :
#6

I expected that when compression is increasing the porosity tends to decrease. But in my case when i increase the Isotropic value - the porosity also increases. In the initial script when you increase pressure the porosity decreases. But when I added the parameter of young module in the O.materials.append it showed the wrong tendency I mentioned before

Sent from my iPhone

> On 24 Nov 2021, at 00:20, Jan Stránský <email address hidden> wrote:
>
> Your question #699605 on Yade changed:
> https://answers.launchpad.net/yade/+question/699605
>
> Status: Open => Needs information
>
> Jan Stránský requested more information:
> Thanks for clarification.
>
> What does "inadequate" mean here in the relation of porosity vs. young?
> What is expected behavior?
>
> Cheers
> Jan
>
> --
> To answer this request for more information, you can either reply to
> this email or enter your reply at the following page:
> https://answers.launchpad.net/yade/+question/699605
>
> You received this question notification because you asked the question.

Revision history for this message
Jan Stránský (honzik) said :
#7

> porrr=utils.voxelPorosity(200,(0.75,0.75,0.75),(1.25,1.25,1.25))

The problem here is voxelPorosity and its parameters.
resolution=200 might be too little, you can try higher values.

More importantly, you are using some random region for evaluation, porosity computed there has (or at least may have) no relation to the true total porosity. Randomly omitting / adding some border parts of spheres influence significantly the result.
You can play with shifting the region.

If I use an analytical formula
###
def porosity():
    volTot = O.cell.volume
    r = 0.09
    volSphs = len(O.bodies) * 4/3. * pi * r**3
    ds = [2*r - i.geom.penetrationDepth for i in O.interactions]
    # https://mathworld.wolfram.com/Sphere-SphereIntersection.html
    volContactss = [1/12.*pi*(4*r+d)*(2*r-d)**2 for d in ds]
    volContacts = sum(volContactss)
    volSolid = volSphs - volContacts
    volVoid = volTot - volSolid
    return volVoid / volTot
###
I get consistent results:
0.47577 (isotropic pressure -1e5)
0.47326 (isotropic pressure -5e5)
0.47202 (isotropic pressure -7e5)

Cheers
Jan

Revision history for this message
Fedor (vodovozov) said :
#8

Thank you very much , but I thought that orthogonal packing ( max porosity) would have turned into hexagonal packing ( min porosity) during compression and the porosity would decrease. But now I realized that change of porosity occurs thank to intersection (overlapping) particles. Does my suggestion make sense?
But in the real soil porosity changes due to particles compaction not intersection.
But again when I started working with initial script I obtained good relation pressure/porosity but when I added E , I was confused.

Sent from my iPhone

> On 24 Nov 2021, at 01:41, Jan Stránský <email address hidden> wrote:
>
> Your question #699605 on Yade changed:
> https://answers.launchpad.net/yade/+question/699605
>
> Status: Open => Answered
>
> Jan Stránský proposed the following answer:
>> porrr=utils.voxelPorosity(200,(0.75,0.75,0.75),(1.25,1.25,1.25))
>
> The problem here is voxelPorosity and its parameters.
> resolution=200 might be too little, you can try higher values.
>
> More importantly, you are using some random region for evaluation, porosity computed there has (or at least may have) no relation to the true total porosity. Randomly omitting / adding some border parts of spheres influence significantly the result.
> You can play with shifting the region.
>
> If I use an analytical formula
> ###
> def porosity():
> volTot = O.cell.volume
> r = 0.09
> volSphs = len(O.bodies) * 4/3. * pi * r**3
> ds = [2*r - i.geom.penetrationDepth for i in O.interactions]
> # https://mathworld.wolfram.com/Sphere-SphereIntersection.html
> volContactss = [1/12.*pi*(4*r+d)*(2*r-d)**2 for d in ds]
> volContacts = sum(volContactss)
> volSolid = volSphs - volContacts
> volVoid = volTot - volSolid
> return volVoid / volTot
> ###
> I get consistent results:
> 0.47577 (isotropic pressure -1e5)
> 0.47326 (isotropic pressure -5e5)
> 0.47202 (isotropic pressure -7e5)
>
> Cheers
> Jan
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/699605/+confirm?answer_id=6
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/699605
>
> You received this question notification because you asked the question.

Revision history for this message
Jan Stránský (honzik) said :
#9

> but I thought that orthogonal packing would have turned into hexagonal packing

it remains orthogonal.
Under the loading conditions, the packing seems to be stable like that.

> But now I realized that change of porosity occurs thank to intersection (overlapping) particles. Does my suggestion make sense?

yes.
(There have been some discussions about interpretation of such results and that porosity computed in this way may not have physical meaning, but it is another topic)

> But in the real soil porosity changes due to particles compaction not intersection.

yes, a more realistic packing should behave like that

> But again when I started working with initial script I obtained good relation pressure/porosity

what is "initial script"? The one you linked?
If so, there is hexagonal packing. As the periodic cell is cubic, the space is not fully filled and the spheres have some freedom to move, unlike in the case of the orthogonal packing.

> but when I added E , I was confused.

E is of course present (although it is hidden) also in the original script, with the default value 1e7.

Cheers
Jan

Revision history for this message
Fedor (vodovozov) said :
#10

Thank you very much, but whether I use makeCloud packing, can I avoid the particles intersections or not?
Big thanks once again for your previous elaborate reply

Sent from my iPhone

> On 24 Nov 2021, at 02:25, Jan Stránský <email address hidden> wrote:
>
> Your question #699605 on Yade changed:
> https://answers.launchpad.net/yade/+question/699605
>
> Status: Open => Answered
>
> Jan Stránský proposed the following answer:
>> but I thought that orthogonal packing would have turned into hexagonal
> packing
>
> it remains orthogonal.
> Under the loading conditions, the packing seems to be stable like that.
>
>> But now I realized that change of porosity occurs thank to
> intersection (overlapping) particles. Does my suggestion make sense?
>
> yes.
> (There have been some discussions about interpretation of such results and that porosity computed in this way may not have physical meaning, but it is another topic)
>
>> But in the real soil porosity changes due to particles compaction not
> intersection.
>
> yes, a more realistic packing should behave like that
>
>> But again when I started working with initial script I obtained good
> relation pressure/porosity
>
> what is "initial script"? The one you linked?
> If so, there is hexagonal packing. As the periodic cell is cubic, the space is not fully filled and the spheres have some freedom to move, unlike in the case of the orthogonal packing.
>
>> but when I added E , I was confused.
>
> E is of course present (although it is hidden) also in the original
> script, with the default value 1e7.
>
> Cheers
> Jan
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/699605/+confirm?answer_id=8
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/699605
>
> You received this question notification because you asked the question.

Revision history for this message
Best Jan Stránský (honzik) said :
#11

> but whether I use makeCloud packing, can I avoid the particles intersections or not?

intersections are essential part of DEM, you cannot avoid it..

If you can omit it in your porosity computations, yes for "normal" circumstances.
As you see from your original results and analytical results, the error of voxelPorosity may be of orders of magnitude larger then the intersections.

Cheers
Jan

Revision history for this message
Fedor (vodovozov) said :
#12

Thank you for your clarifications.

Revision history for this message
Fedor (vodovozov) said :
#13

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