# orthogonal packing

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:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Robert Caulk (rcaulk) said on 2021-11-23: #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

 Revision history for this message Fedor (vodovozov) said on 2021-11-23: #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
PeriTriaxController(
label='triax',
# specify target values and whether they are strains or stresses
goal=(sigmaIso, sigmaIso, sigmaIso),
# 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
# 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)
# 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 on 2021-11-23: #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 on 2021-11-23: #4

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:
>
> Status: Open => Needs 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
>
> --
> this email or enter your reply at the following page:
>
> You received this question notification because you asked the question.
>

 Revision history for this message Jan Stránský (honzik) said on 2021-11-23: #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 on 2021-11-23: #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:
>
> Status: Open => Needs information
>
> Thanks for clarification.
>
> What does "inadequate" mean here in the relation of porosity vs. young?
> What is expected behavior?
>
> Cheers
> Jan
>
> --
> this email or enter your reply at the following page:
>
> You received this question notification because you asked the question.

 Revision history for this message Jan Stránský (honzik) said on 2021-11-23: #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 on 2021-11-23: #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:
>
> 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:
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
>
> You received this question notification because you asked the question.

 Revision history for this message Jan Stránský (honzik) said on 2021-11-23: #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 on 2021-11-23: #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:
>
> 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:
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
>
> You received this question notification because you asked the question.

 Revision history for this message Jan Stránský (honzik) said on 2021-11-23: #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 on 2021-11-24: #12

Thank you for your clarifications.

 Revision history for this message Fedor (vodovozov) said on 2021-11-24: #13

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

To post a message you must log in.