Yade

Cannot simulate periodic box with walls

Asked by Danny (Dingeman) van der Haven

Dear YADE users,

My goal is to simulate a representative volume element under uniaxial compression. I want to compare two cases:
1) a fully periodic simulation in which the particle coordinates are scaled
2) a simulation in which the particles are compressed by walls at the top and bottom of the box, with periodicity in all other directions.

Case 1 is easy to simulate, using tutorial example 06-periodic-triaxial-test.py. To simulate case 2 I used example 06-periodic-triaxial-test.py as a starting point and started integrating elements from 03-oedometric-test.py.

I can insert a wall in the periodic box of example 06. However, as soon as I also add the sphere-wall interaction to the collider, the simulation no longer runs. YADE gives no error messages but refuses to iterate past iteration 1.

The input script I would like to run is (note that although this script is more complicated, I already run into issues when only adding the wall and sphere-wall interactions to example 06):

from __future__ import print_function
sigmaIso = -1e5
sigmaMin = -1e3

#import matplotlib
#matplotlib.use('Agg')

# generate loose packing
from yade import pack, qt, plot

O.periodic = True
sp = pack.SpherePack()
# uniform distribution
sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.3, periodic=True)

# setup periodic boundary, insert the packing
sp.toSimulation()

O.bodies.append(wall(0.0, axis=2, sense=1))

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()], allowBiggerThanPeriod=True),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
                        [Ip2_FrictMat_FrictMat_FrictPhys()],
                        [Law2_ScGeom_FrictPhys_CundallStrack()]),
        NewtonIntegrator(damping=.2),
        PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker')
]
O.dt = .25 * PWaveTimeStep()

def checkUnbalanced():
        if O.iter < 5000:
                return
        if unbalancedForce() > 0.1:
                return
        O.bodies.append(wall(2.0, axis=2,sense=-1))
        global plate
        plate.state.vel = (0,0,-0.01)
        O.engines = O.engines + [PyRunner(command='addPlotData()',iterPeriod=200)]
        checker.command = 'unloadPlate()'

def unloadPlate():
        if abs(O.forces.f(plate.id)[2]/(2*2)) > abs(sigmaIso):
                plate.state.vel *= -1
                checker.command = 'stopUnloading()'

def stopUnloading():
        if abs(O.forces.f(plate.id)[2]/(2*2)) < abs(sigmaMin):
                plot.saveDataTxt(O.tags['d.id']+'.txt')
                O.pause()

def addPlotData():
        if not isinstance(O.bodies[-1].shape,Wall):
                plot.addData()
                return
        Fz = O.forces.f(plate.id)[2]
        plot.addData(
                unbalanced=unbalancedForce(),
                i=O.iter,
                sz=abs(fz)/(2*2),
                ev=abs(plate.state.pos[2]*2*2-2*2*2)/(2*2*2)
        )

# define what to plot
plot.plots = {'i': ('unbalanced',),'ev': ('sz'),}
# show the plot
plot.plot()

O.run()

With kind regards,
Danny van der Haven

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
Jan Stránský (honzik) said : 		#1

Hello,

> 2) a simulation in which the particles are compressed by walls at the top and bottom of the box, with periodicity in all other directions.

Why not simply let periodicity to control also the "top-bottom" direction? Using PeriTriaxController or similar existing tools. Are walls necessary?

There is obviously some problem with walls and this simulation settings [1]. Have a try with a box, like at [2]

Cheers
Jan

[1] https://answers.launchpad.net/yade/+question/692092
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/PeriodicBoundaries/periodicSandPile.py

Revision history for this message
Danny (Dingeman) van der Haven (dlhvdh) said : 		#2

Dear Jan,

Thank you for your answer.

Some sort of wall is essential. Correct me if I am wrong, but when the box is changed using PeriTriaxController it seems that all particle coordinates are scaled. This means that stresses will be distributed uniformly. This is not the case in real uniaxial compaction, where the stress varies depending on the height of the slice you are looking at. For this reason, I want to simulate a case where the pressure is increased using plates.

Using a box as given in [2] did not work, PeriTriaxController does not scale the inserted box geometry. A hacky alternative that seems to work is to create two thin boxes (one at each end) and use those as plates. However, if I fix these plates and use PeriTriaxController the particle coordinates are still scaled. If I do not fix the plates and give them a velocity then the plates move away. How can I make the box-plates move at a constant speed (irrespective of the other forces on them)?

Also, the documentation says the option dynCell determines how the imposed stress is controlled. What is meant with "packing stiffness" and "applying the laws of dynamics"? I assume the latter computers the average forces on the particles and uses that to compute a stress.

Please let me know if I should add a bug report on GitLab about the problem with walls and periodic cells. Similarly, if there is a way that I can help improve the documentation on the YADE website, please let me know. I plan on implementing new features in YADE and would like to get involved with the development.

With kind regards,
Danny

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

> Some sort of wall is essential

depends on what "essential" is.
Also depends on what "some sort" is, because periodic boundary conditions essentially are "some sort of wall"..

> when the box is changed using PeriTriaxController it seems that all particle coordinates are scaled.

basically yes, but it can be controlled and changed, see [3].
If you set homoDeform=0, you should not get the "scaling"

> This means that stresses will be distributed uniformly.

rather strain is distributed uniformly, stress is then result of stiffness and packing structure.
If the packing is "uniform", then also the resulting stress would be "uniform".
However, if the packing is "non-uniform", also the resulting stress would be "non-uniform".

The "scaling" is only one component of the particles motion, and therefore only one component to the resulting stress. There are of course other components (interparticle interactions mainly) influencing the result.

> This is not the case in real uniaxial compaction, where the stress varies depending on the height of the slice you are looking at.

What is the reason of stress variation? Dynamics? Gravity? Something else?
For "fast" loading yes, the stress waves should not be neglected and "wall" loading is more realistic.
Contrary, for "slow" loading, this "periodic approach" is actually more realistic than the loading induced only by planes (which induce stress waves not present in reality).

> Using a box as given in [2] did not work

is it still freezing, or just the behavior and results are not as expected?

> PeriTriaxController ...

I suggest either:
- use PeriTriaxController with no walls
or:
- use "enlarged" periodic cell with only two periodic (transverse) directions and the axial direction enlarged (essentially making it non-periodic axially). Then you can push the sample by walls and somehow control stress/strain in the transverse directions.
or:
- probably both approaches can be used together. The plates are moving "independently", the axial direction is strain-controlled with zero change, the transverse directions are controlled by PeriTriaxController. Have a try :-)

> If I do not fix the plates and give them a velocity then the plates move away.
> How can I make the box-plates move at a constant speed (irrespective of the other forces on them)?

b.state.blockedDOFs = "xyzXYZ # [4]
b.state.vel = someConstantVelocity

> ... dynCell

basically yes, you can compute stress assuming static equilibrium or you can somehow take into account inertia of the particles.

> Please let me know if I should add a bug report on GitLab about the problem with walls and periodic cells.

yes, it seems like a bug

> if there is a way that I can help improve the documentation on the YADE website, please let me know.

Yes. the documentation is part of the source code. You "just" download the source code, modify it, and commit the changes to the gitlab. I don't remember exactly the process, but it is described in documentation and is out of scope of this question anyway (of course you are free to open a new question aiming on this).

I plan on implementing new features in YADE and would like to get involved with the development.

Cool, good luck :-)

Cheers
Jan

[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Cell.homoDeform
[4] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.State.blockedDOFs

Revision history for this message
Danny (Dingeman) van der Haven (dlhvdh) said (last edit ): 		#4

Dear Jan,

> depends on what "essential" is.
> Also depends on what "some sort" is, because periodic boundary conditions essentially are "some > sort of wall"..

A wall that interacts with the particles like a physical wall is necessary for me to be able to make the comparison I mentioned in the original question.

> What is the reason of stress variation? Dynamics? Gravity? Something else?
> For "fast" loading yes, the stress waves should not be neglected and "wall" loading is more realistic.
> Contrary, for "slow" loading, this "periodic approach" is actually more realistic than the loading
> induced only by planes (which induce stress waves not present in reality).

These are valid concerns / questions. This is why I want to do the comparison between plate compression and periodic compression. I study the compaction of pharmaceutical powders. The speed of the compression can be 0.1 mm/s (quasi-static) but also 950.00 mm/s (dynamic) for tablets with a final thickness of 4.00 mm. Gravity is not a concern, but dynamics and wall interactions are. For slow compression speeds the variation in stress along the compression axis is likely caused by interaction with the side walls. The comparison is mainly to look at the possible differences made by dynamics.

> is it still freezing, or just the behavior and results are not as expected?

The simulation iterates, but the behaviour is not what I was aiming for. The box stays the same size so the wall do not exert any pressure on the particles. The cell does change.

> If you set homoDeform=0, you should not get the "scaling"

Thank you, that option indeed disabled the scaling / shifting of particle coordinates.

> b.state.blockedDOFs = "xyzXYZ”
> b.state.vel = someConstantVel

This solved my problem! Thank you! I created a periodic box with an extra long z-axis. Using two box-plates I can then compress the particles in the z-direction using the blockedDOFs and vel options, keeping periodicity in the x and y directions.

> basically yes, you can compute stress assuming static equilibrium or you can somehow take into account inertia of the particles.

So that means if dynCell = False you assume static equilibrium and if dynCell = True you take into account the inertia?

> yes, it seems like a bug

I filed an issue report on Gitlab.

Kind regards,
Danny

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

concerning dynCell, I don't use it myself and don't know much details.
If you are interested on more details, please open a new question focusing on dynCell
Cheers
Jan

Revision history for this message
Jérôme Duriez (jduriez) said : 		#6

Regarding dynCell, maybe Eq. (47) of https://www.sciencedirect.com/science/article/pii/S1631072117302346 could help to understand what are these "laws of dynamic" here.

(Including this reference is actually on my TODO list for Periodic Boundaries documentation, following up in some future on https://gitlab.com/yade-dev/trunk/-/merge_requests/802)

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#7

Hi Danny,
I didn't read all posts in details but I suspect Jan already suggested how to proceed: insert boxes and impose their velocities.
Meanwhile the periodic cell should not move.
If the problem persist please post your script again.
Cheers
Bruno

ps. I closed the gitlab issue [2] since it is not actually a bug in my view (or at least, it won't be fixed).

[2] https://gitlab.com/yade-dev/trunk/-/issues/246

Revision history for this message
Danny (Dingeman) van der Haven (dlhvdh) said : 		#8