quasi-static equilibrium

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://answers.launchpad.net/yade/+question/136034
>
> 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://answers.launchpad.net/yade/+question/136034
>
> 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://answers.launchpad.net/yade/+question/136034
> >
> > 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://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
> More help : https://help.launchpad.net/ListHelp
>

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={Rheophysics of dense granular materials{:} discrete simulation
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://answers.launchpad.net/yade/+question/136034
>
> 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

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.unabalncedForce()>0.05), which is unavoidable, seems to me. The subsequent evolution, though, as interactions get established, could be made more corresponding to expectations by multiplying e.g. by Mi/M (where Mi would be mass of particles which are in interactions, and M total mass of all particles) or something similar.

(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/control is one that will _not_ affect the
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>/msg00071.html.
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.

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.

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-theorem-like form: energy flux over the boundary (influx) = energy dissipation (efflux, or divergence) + free energy in the system. Note that we know the influx (boundary control), and if we can estimate the dissipation term, then the free energy estimate is given "for free" (kinetic energy, plastic potential etc etc).

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→dynamic→static, where in the first time the influx is bigger than dissipationv in the initial "dynamization", and it inverses for the final stabilization. The question is how much "dynamic" we want to get in the middle phase -- it could be given by some energy measure, such as maximum admissible energy concentration [Jm⁻³], or enegy concentration per interaction area [Jm⁻³m⁻²].

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-than-necessary dynamization, and slower-than-necessary stabilization).

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 dynamization, and 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

> 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