# quasi-static equilibrium

Hello,

I would like to address you a question about static-equilibrium in dem simulations, with particular reference to PeriTriaxController engine. In the latter, how do we make sure that the static-equilibrium condition is maintained during the triaxial phase, for instance? Looking at the code, I see that the unbalanced force is checked only when the system has reached the prescribed goal condition. What is the best criterion to establish if the system is under quasi-static conditions? Is the unbalanced force the only solution?

Thank you for your opinion.

cheers, Chiara

## Question information

- Language:
- English Edit question

- Status:
- Answered

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Last query:
- 2010-11-30

- Last reply:
- 2010-12-01

Hi Chiara,

Not the only one, but unbalanced force is the best solution to my

knowledge, since it really evaluate out-of-staticity forces in a

non-dimensional form.

I always used that in papers and never got a comment.

I've seen disadvantages in all other indicator proposed (except variants

of unbalanced force, like fb_max/fc_mean - you can get this one with a

flag).

Bruno

Le 30/11/2010 10:40, Chiara Modenese a écrit :

> New question #136034 on Yade:

> https:/

>

> Hello,

>

> I would like to address you a question about static-equilibrium in dem simulations, with particular reference to PeriTriaxController engine. In the latter, how do we make sure that the static-equilibrium condition is maintained during the triaxial phase, for instance? Looking at the code, I see that the unbalanced force is checked only when the system has reached the prescribed goal condition. What is the best criterion to establish if the system is under quasi-static conditions? Is the unbalanced force the only solution?

>

> Thank you for your opinion.

> cheers, Chiara

>

>

>

>

> Ce message entrant est certifie sans virus connu.

> Analyse effectuee par AVG - www.avg.fr

> Version: 9.0.872 / Base de donnees virale: 271.1.1/3288 - Date: 11/29/10 20:34:00

>

Hi all! Thanks for your answers! I see that the unbalanced force is the best

indicator. Another matter is how to calibrate the algorithm's parameters to

achieve quasi-static equilibrium. I am used to scale up the mass, but when I

employ high inter-particle stiffnesses, then the system is quite unstable

and reducing the time step does not help. There are so many parameters

involved, like mass, loading rate, damping and so on. I could not find a

good paper where they suggest a "precise" way of calibration to reach

quasi-static equilibrium. For instance, can we use a smaller sample to

calibrate those parameters? Would be the calibration still representative?

cheers! Chiara

On 30 November 2010 11:30, Chareyre <email address hidden>wrote:

> Question #136034 on Yade changed:

> https:/

>

> Status: Open => Answered

>

> Chareyre proposed the following answer:

> Hi Chiara,

>

> Not the only one, but unbalanced force is the best solution to my

> knowledge, since it really evaluate out-of-staticity forces in a

> non-dimensional form.

> I always used that in papers and never got a comment.

> I've seen disadvantages in all other indicator proposed (except variants

> of unbalanced force, like fb_max/fc_mean - you can get this one with a

> flag).

>

> Bruno

>

> Le 30/11/2010 10:40, Chiara Modenese a écrit :

> > New question #136034 on Yade:

> > https:/

> >

> > Hello,

> >

> > I would like to address you a question about static-equilibrium in dem

> simulations, with particular reference to PeriTriaxController engine. In the

> latter, how do we make sure that the static-equilibrium condition is

> maintained during the triaxial phase, for instance? Looking at the code, I

> see that the unbalanced force is checked only when the system has reached

> the prescribed goal condition. What is the best criterion to establish if

> the system is under quasi-static conditions? Is the unbalanced force the

> only solution?

> >

> > Thank you for your opinion.

> > cheers, Chiara

> >

> >

> >

> >

> > Ce message entrant est certifie sans virus connu.

> > Analyse effectuee par AVG - www.avg.fr

> > Version: 9.0.872 / Base de donnees virale: 271.1.1/3288 - Date: 11/29/10

> 20:34:00

> >

>

> --

> You received this question notification because you are a member of

> yade-users, which is an answer contact for Yade.

>

> _______

> Mailing list: https:/

> Post to : <email address hidden>

> Unsubscribe : https:/

> More help : https:/

>

Luc Sibille (luc-sibille) said : | #3 |

Sorry, I have forgotten tne new way to reply to questions. Thank you Anton.

Luc

Hi Chiara,

In a general way, I think that one of the best answers is to read to PhD

Thesis of Bruno where he checked in which conditions we can assume that

the solution given by the DEM (which is a dynamic solution) for triaxial

compressions is not very far from a static solution. However his thesis

is written in French ;-)

Roughly speaking, the idea is to show that for sufficiently small

loading velocities (but not too small otherwise the simulation takes to

much time), the solution given by the DEM is lowly influenced by the

loading velocity. The "lowly" is fixed by what you think to be

sufficiently accurate for the problem you study.

Then for practical and repetitive simulations you can check the

unbalance force and verify that it corresponds to the unbalance force

obtained during the parametric study to find the "sufficiently low

loading velocity".

Maybe Bruno will tell you that you can directly reason from the

unbalance force?

A slightly different approach, but consisting also to estimate the

importance of the inertial terms through the dimensionless strain rate I

can be found in:

@ARTICLE{dacruz05,

author={da Cruz, F. and Emam, S. and Prochnow, M. and Roux, J.N. and

Chevoir, F.},

title=

of plane shear flows},

journal={Physical Review E},

volume={72},

number={2},

pages = "021309",

year={2005},

note={},

}

Chiara Modenese a écrit :

> New question #136034 on Yade:

> https:/

>

> Hello,

>

> I would like to address you a question about static-equilibrium in dem simulations, with particular reference to PeriTriaxController engine. In the latter, how do we make sure that the static-equilibrium condition is maintained during the triaxial phase, for instance? Looking at the code, I see that the unbalanced force is checked only when the system has reached the prescribed goal condition. What is the best criterion to establish if the system is under quasi-static conditions? Is the unbalanced force the only solution?

>

> Thank you for your opinion.

> cheers, Chiara

>

--

Luc Sibille

Université de Nantes - Laboratoire GeM UMR CNRS

IUT de Saint Nazaire

58, rue Michel-Ange - BP 420

44606 Saint-Nazaire Cedex, France

Tel: +33 (0)2 40 17 81 78

Václav Šmilauer (eudoxos) said : | #4 |

I think the criterion be derived from an exact definition of "quasi-static", which we don't have (the definition is necessarily also "quasi"). It is clear that it is a measure that will be zero for static case (no motion, or more generally no accelerations, since we are in newtonian physics); since exact zero will be only rarely achieved, one needs to set some tolerance, and therefore it should be an admimensional measure, since otherwise the threshold will have different meanings for different problem scale.

The easiest adimensional params are ratios of 2 quantities with the same dimension; unbalanced force is one, but there are other possible: I am not sure if energy balance would make any difference (kinetic/elastic potential) -- it could, since physical systems conserve energy rather than forces, therefore even if they are both adimensional and one will be monotnically increasing function of the other, it can be a better measure. Never tried, though.

The disadvantage of unbalanced force (or the energy ratio) is that it is not defined when there are no interactions (leading to conjuncted conditions like O.iter>1000 && utils.unabalnce

(I am sorry to have no good answer, take it just as a comment...)

> not sure if energy balance would make any difference (kinetic/elastic

> potential) -- it could, since physical systems conserve energy rather

> than forces, therefore even if they are both adimensional and one will

> be monotnically increasing function of the other, it can be a better

> measure. Never tried, though.

This is another reasonable option. The problem is you can have free particles in holes of

the packing, spining forever since they have no contact. They will give a finite kinetic

energy, but null unbalanced force. The latest is more relevant since such particles do not

correspond to non-static behaviour of the packing.

Bruno

> achieve quasi-static equilibrium. I am used to scale up the mass, but when I

> employ high inter-particle stiffnesses, then the system is quite unstable

> and reducing the time step does not help. There are so many parameters

> involved, like mass, loading rate, damping and so on.

1- Is the system stable (i.e. unbalanced force monotonicaly decreasing down to 1e-10 or

so) under isotropic (or gravitational - I don't know what you are doing exactly) loading?

2- If yes, do you find a proportionality between loading rate and unbalanced force during

the loading?

3- If yes, unbF/rate is most probably scaled with damping. Do you really need damping?

4- If yes, find the optimal damping in terms of CPU time (most probably around 0.1-0.2).

5- else get rif of damping

6- else, there is a problem in the formulation (either contact law or boundary

conditions control). It needs to be fixed.

7- else, same as 6.

> I could not find a

> good paper where they suggest a "precise" way of calibration to reach

> quasi-static equilibrium. For instance, can we use a smaller sample to

> calibrate those parameters? Would be the calibration still representative?

I don't think such paper exists, especially if you use Hertz-Mindlin.

Generally, no, you will not get exactly the same behaviour with a different number of

particles. Numerical waves have a velocity of 1 particle/timestep in explicit scheme,

whatever dt and size of particle. Hence they will be different with different size of packing.

OTOH, if there is a problem in the formulation (6,7), it is a good idea to investigate it

in small packings.

In triaxial tests, I found different ratio between unbalanced force and peak stress error

with different number of particles.

Bruno

There was an indentation bug, sorry. ;-)

It should read :

1- Is the system stable (i.e. unbalanced force monotonicaly decreasing down to 1e-10 or

so) under isotropic (or gravitational - I don't know what you are doing exactly) loading?

2- If yes, do you find a proportionality between loading rate and unbalanced force during

the loading?

3- If yes, unbF/rate is most probably scaled with damping. Do you really need damping?

4- If yes, find the optimal damping in terms of CPU time (most probably around 0.1-0.2).

5- else get rif of damping

6- else, there is a problem in the formulation (either contact law or boundary

conditions control). It needs to be fixed.

7- else, same as 6.

Bruno

Thanks for your answers. I do not expect exactly the same solution for a reduced number of particles, but maybe the parameters calibrated for it would still work once the size of the simulation is increased.

To my knowledge, no paper on this sort of calibration for a non-linear elastic law (usually authors say that a trial and error procedure is adopted, but not so much details on it are provided - OTOH, if one wants to reproduce the same results, such info are quite important to prove the correctness of those results, I think).

My sample gets the stability after the isotropic phase. The instability arises when the triaxial path starts and only when I increase the contact Young's modulus. I cannot get rid of damping. I was looking for a new procedure to control the boundaries in order to solve this problem. Not sure yet what to adopt. I do not think that the current procedure (the one with the dynamic cell) is bad, but I would like to find a more systematic manner to approach the problem.

Chiara

> OTOH, if one wants to reproduce the same results, such info are quite important to prove the correctness of those results, I think).

Not really. A good setting of rate/damping/

result. Not affecting means a small change in the settings keeps the results statisticaly

unchanged.

Hence, if Thornton used settings A and you use settings B, and if nor A nor B affect the

result. It is consistent to compare the results.

> My sample gets the stability after the isotropic phase. The instability arises when the triaxial path starts and only when I increase the contact Young's modulus. I cannot get rid of damping. I was looking for a new procedure to control the boundaries in order to solve this problem. Not sure yet what to adopt.

If boundary control is stable for isotropic confinment, it is unlikely that it is the

source of instabilities in the loading phase (provided the loading rate is low enough).

Do you increase Young during the simulation? It could explain a lot.

Bruno

I meant: I would like to see a work in which authors also present the calibration analysis in itself (for the non-linear elastic contact case, for instance), where the influence of the parameters (and their meanings) involved is accessed. Then, of course, different procedures lead to the same results as long as they do not affect them.

I keep the Young's modulus the same during simulation. When I change its value I refer to another sample. Actually I realize now that perhaps I am not setting the max unbalanced force low enough :-( I see you mentioned in this post a value for it of 1e-10. so maybe my system is, in fact, dynamic. The I just found this interesting discussion you had some time ago: http://<email address hidden>

I think I need more analyses to see what happens.

Anyway thanks for your suggestions so far. Cheers! Chiara

In my case, am also using contact viscous damping. Though this should not have an influence, I think, since the condition is quasi-static. But it should help to reach equilibrium. But again, I have not tried yet to vary it, so am not sure how it really works.

Chiara

> I meant: I would like to see a work in which authors also present the calibration analysis in itself (for the non-linear elastic contact >case, for instance), where the influence of the parameters (and their meanings) involved is accessed. Then, of course, different >procedures lead to the same results as long as they do not affect them.

Then read my thesis :-) (thanks Luc for kind aknowledgement!)

> I keep the Young's modulus the same during simulation. When I change its value I refer to another sample. Actually I realize now that >perhaps I am not setting the max unbalanced force low enough :-( I see you mentioned in this post a value for it of 1e-10.

No, don't consider 10e-10 a target value in routine simulations, or you will wait forever.

It is only a good test to check that unbF is monotonically decreasing with steps in the

asbence of loading fluctuations (and past 10e-10, don't expect further decrease because

you reached numerical noise). 10e-3 is a reasonable default threshold in TT in my opinion

(but I'm quite liberal on this aspect ;-)).

Cheers.

B.

Jérôme Duriez (jduriez) said : | #13 |

Le 30/11/2010 15:39, Chareyre a écrit :

> 2- If yes, do you find a proportionality between loading rate and unbalanced force during

> the loading?

>

>

This is new for me. From where does it come from ? How do you compare

loading rate (strain rate for example) with unbalanced force ?

> This is new for me. From where does it come from ? How do you compare

> loading rate (strain rate for example) with unbalanced force ?

Multiply the loading rate by 2 and see the effect on average unbF in the simulation?

Of course, I'm not speaking of exact proportionality here (though it was quite linear in

the range I tested).

B.

Václav Šmilauer (eudoxos) said : | #15 |

A few thoughts on this one, after some discussion:

1. PeriTriaxController confounds (in the mass param) static and dynamic response of the packing, whereas they could be separated in principle. By static, I mean stiffness, i.e. what stress delta corresponds to what strain delta (for stress-control); to my knowledge, this is currently computed well by Peri3dController (surrent global tangent stiffness).

2. The measure of staticity depends on the formulation used (damping, for instance); the example Bruno gave (with free sphere rotating) is perhaps good for cundall damping, though I would personally try to employ another daming algorithm to get rid of the energy anyway (imagine the sphere touches at some point another particle, and the energy is back there -- increase in unbalanced force witout any work on the boundary).

Intergral (over the whole system) energy measures are consistent, since energy is what is conserved in the system. If one can formulate how much energy can be dissipated at each simulation stage, then this can be compared to boundary work. Currently, we compare apples and oranges (maxStrainRate and mass for work, numerical damping for dissipation).

3. The energy "balance" between work on the boundary (increment of meanfield kinetic energy approx, which immediately transforms in elastic potential and fluctuation kinetic energy; real work on boundary for the aperiodic case) and the dissipation (plasticity, damping, damage, ...) can be written in a divergence-

4. If the process starts from static state (supposition), then it cannot be static all the time (nothing would move); we need to describe the transition static→

The transition could be controlled by some smooth function f(t=0)=0, f(t=1)=1, where t=0 is the intial static state and t=1 is the final static state. Then f(t) would give relative predominance of dissipation over the influx (and vice versa for 1-f(t) obviously).

Currently, only maxStrainRate controls the proportion; it is dimensionally incorrect (wherefore the necessity for trial-and-error) and constant (it is likely to give too low influx at the beginning, and too high at the end -- slower-

A few comments :

1.a) PeriTriaxController (with dynCell=True) is not confounding statics and dynamics : it

is purely dynamic.

1.b) I recall there is no exact theoretical derivation of stiffness from contact list for

"granular" materials (i.e. 3-4 interactions per element, no cohesion nor rolling

friction); and that is the reason why I didn't use for periodic case the static strategy

used for TT with boxes.

> 2. The measure of staticity depends on the formulation used (damping,

> for instance);

I don't see why unbF would depend on formulation, since it only employs forces.

> the example Bruno gave (with free sphere rotating) is

> perhaps good for cundall damping, though I would personally try to

> employ another daming algorithm to get rid of the energy anyway (imagine

> the sphere touches at some point another particle, and the energy is

> back there -- increase in unbalanced force witout any work on the

> boundary).

When the sphere touch back, you will have a small increase in unbF. It seems consistent :

you lost a bit of staticity. I see no problem in that. It is not really different than

breaking a contact at iteration N, so that unbF(N+1)>unbF(N) due to ongoing internal

reorganisations.

> Intergral (over the whole system) energy measures are consistent, since

> energy is what is conserved in the system.

Energy measure is consistent, but evaluating staticity and energy conservation are two

different things. Kinetic vs. elastic energy ratio is a good measure of staticity, but I

still have to hear what is the disadvantage of unbF, I don't see any.

Considering that tracing energy is cpu-time consuming, and that there are so many laws,

damping methods, and boundary conditions for which energy tracing is not implemented, and

will not any time soon due to excessive complexity, I don't see energy as a candidate for

routine staticity evaluation.

OTOH, unbF applies for any model where forces are defined, i.e. everything, without

specific lines of code.

> Currently, we compare apples and oranges

> (maxStrainRate and mass for work, numerical damping for dissipation).

I don't understand what comparisons are refered to here. Out-of-staticity errors and

damping (or dissipation) are two different things, even if they interact on each other.

You can have out-of-staticity effects, without damping, in a purely elastic problem :

apply a force on a steel bar at t=0, get the displacement to define the force/disp

relation for the bar. Ouch : disp. is a sinusoidal function of time!

> 3. The energy "balance" between work on the boundary (increment of

> meanfield kinetic energy approx, which immediately transforms in elastic

> potential and fluctuation kinetic energy; [...]

This is all good, but I don't think you are really speaking of non-staticity errors here.

Whatever you do with energy, you will not escape the problem of numerical waves : a

perturbation (typically contact breakage inside the packing) trigger non-equilibrium. It

will need a number of steps before the system goes back to equilibrium. This number

depends on physics, but also on the explicit FD scheme that propagates information at the

speed of 1 particle/iteration.

In fact, in a triaxial simulation of ~20k spheres, such perturbations are constantly

appearing and stabilizing (making sure that they all occure independently one after the

other would need weeks-long simulations). You have to choose a loading rate so that

stabilization process will be fast enough compared to loading, and this is not a matter of

energy.

I agree with what you say about energy fluxes and all. That is really interesting, but not

practical for choosing a loading rate.

Most of the quasi-static stabilisation processes have a velocity in particle/iteration,

not in m/s, that is the key point (same for loading rate, that can be considered in s-1 or

iter-1 depending on what effects matter).

> Currently, only maxStrainRate controls the proportion; it is

> dimensionally incorrect (wherefore the necessity for trial-and-error)

> and constant (it is likely to give too low influx at the beginning, and

> too high at the end -- slower-

> than-necessary stabilization).

I used to define loading rates with ramps. I stopped bothering with that : in a typical

triaxial loading, the max rate will be imposed by the error you admit on peak and residual

stresses, i.e. typically 0.01<eps<0.3. Increasing rate in the very early loading - if

admissible - will make you save time on a very short interval : around 3% of the total

simulation. It is not worth the effort to calibrate separately the different stages.

Bruno

Václav Šmilauer (eudoxos) said : | #17 |

> 1.a) PeriTriaxController (with dynCell=True) is not confounding statics and dynamics : it

> is purely dynamic.

Eh, that is exactly what I said... It has no notion of static response,

counfindg the two together (if you don't agree, it is a matter of words,

I think)

> 1.b) I recall there is no exact theoretical derivation of stiffness from contact list for

> "granular" materials (i.e. 3-4 interactions per element, no cohesion nor rolling

> friction); and that is the reason why I didn't use for periodic case the static strategy

> used for TT with boxes.

Peri3dController is getting the response from the packing, not from

formulas...

>> 2. The measure of staticity depends on the formulation used (damping,

>> for instance);

> I don't see why unbF would depend on formulation, since it only employs

> forces.

I meant that choice of a suitable staticity criterion depends on what

type of damping is employed.

> I see no problem in that. It is not really different than

> breaking a contact at iteration N, so that unbF(N+1)>unbF(N) due to ongoing internal

> reorganisations.

The point was that free energy is always non-increasing (unless there is

external work), that's why I though it was better.

> Energy measure is consistent, but evaluating staticity and energy conservation are two

> different things. Kinetic vs. elastic energy ratio is a good measure of staticity, but I

> still have to hear what is the disadvantage of unbF, I don't see any.

Above.

> Considering that tracing energy is cpu-time consuming, and that there are so many laws,

> damping methods, and boundary conditions for which energy tracing is not implemented, and

> will not any time soon due to excessive complexity, I don't see energy as a candidate for

> routine staticity evaluation.

Eh, alright.... I think this has become too confused a discussion to be

meaningful.

v.

>> 1.a) PeriTriaxController (with dynCell=True) is not confounding statics and dynamics : it

>> is purely dynamic.

> Eh, that is exactly what I said...

I got it.

>I meant that choice of a suitable staticity criterion depends on what

>type of damping is employed.

Not so, if you mean staticity criterion like "I want error on sigma<10e-3 in the

simulation", and this is how unbF is usefull.

If you set damping higher, or if you modify the damping method (keeping loading rate

constant), it will be reflected in unbF. So, you can adjust the loading rate in order to

keep unbF constant, and expect the same order of error.

Of course, it needs prior tests to decide what unbF correspond to this 10e-3 error on

sigma, but this is required anyway whatever the non-staticity measurment.

> Eh, alright.... I think this has become too confused a discussion to be

> meaningful.

I think the important point is :

Non-staticity is affected by physics but also by numerical artifacts, coming from the

integration scheme (the "numerical wave" part), affecting the relaxation delay, and

affected by the number of grains. The latest is dominant in the quasi-static regime for

large number of grains. Hence, applying energy concepts (or more generally physics) for

this question is not fully adressing the problem, and it will not replace numerical

experiments.

Bruno

## Can you help with this problem?

Provide an answer of your own, or ask Chiara Modenese for more information if necessary.