getStress for ScGeom

Hi Ning,

> New question #199704 on Yade:
> https://answers.launchpad.net/yade/+question/199704
>
> Hi,
>
> I heard you are unifying the stress calculation for YADE. It's really good news. I still have some concern. The getStress() function will return a Matrix3 which is asymmetric. If one want to get symmetric result, he/she can use toVoigt() which returns a Vector6. The conversion between data types is also a little bit annoying. So maybe a flag could be added to the function getStress() to indicate if one wants to symmetrize the stress tensor.

This week we had discussion with Bruno exactly about these topics. The
conclusion is that the asymmetry of periodic cell stress tensor should
be of low order of magnitude (can you check it?). If it is not, it
should be caused by nonzero angular acceleration (and therefore inertial
torques, then the formula used is not valid). Could you please check
also this ( O.forces.t(id) )? One version of new getStress function will
return the symmetric part :-) but if the asymmetry is too large,
something is wrong..

>
> The function I added to stressTensorOfPeriodicCell is calculating the elastic moduli of the packing, the formula is
>
> S_ijkl = 1/V sum_contacts (L^2 (k_n n_i n_j n_k n_l + k_t n_i t_j n_k t_l))
>
> see: (Kuhl et al., 2001) Microplane modelling and particle modelling of cohesive-frictional materials.
> (Luding, 2004) Micro–macro transition for anisotropic, frictional granular packings.
>
> Because both stress and the moduli are summed over all contacts, I combined them to one function to save calculation time.

I was thinking about function estimating overall elastic constants, but
the result depends on constitutive constants of individual links (unlike
in the case of stress), so there should be done some additional work to
make it general function (like to get current elastic stiffness of a
contact in case of nonlinearities etc.)

The time savage is always good :-)

Jan

Hi Jan,

I noticed that you deprecated stressTensorOfPeriodicCell and suggested to use getStress instead. I think this function works for clumps also. Am I right?

To make general expression of elastic moduli, I think we can use equivalent kn and kt in nonlinear contact models (also change L*n_i and L*n_k to L_i and L_k for non-spherical particles).

Ning

On 23 Jun 2012, at 10:26, ceguo wrote:

> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Status: Answered => Open
>
> ceguo is still having a problem:
> Hi Jan,
>
> I noticed that you deprecated stressTensorOfPeriodicCell and suggested
> to use getStress instead. I think this function works for clumps also.
> Am I right?
>
> To make general expression of elastic moduli, I think we can use
> equivalent kn and kt in nonlinear contact models (also change L*n_i and
> L*n_k to L_i and L_k for non-spherical particles).
Hi Ning,
How would you define the equivalent stiffnesses?
Thanks, Chiara

>
> Ning
>
> --
> 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
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Hi Chiara,

Thanks for the reply. The equivalent stiffness is defined as K=dF/dU, where dF is change of contact force and dU is change of overlap or sliding displacement. So e.g. in Hertzian model, equivalent K_n, K_s are (from PFC mannual):

K_n = (2G sqrt(2R)/(3(1-nu)))sqrt(U_n);

K_s = (2(G^2*3*(1-nu)R)^1/3/(2-nu))*F_n^1/3.

So they are not constant but depend on the current overlap and sliding condition.

Ning

OK, then you mean the tangential stiffnesses. I thought you referred to some equivalent constant value.
Thanks for reply,
Chiara

On 23 Jun 2012, at 14:45, ceguo wrote:

> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> ceguo posted a new comment:
> Hi Chiara,
>
> Thanks for the reply. The equivalent stiffness is defined as K=dF/dU,
> where dF is change of contact force and dU is change of overlap or
> sliding displacement. So e.g. in Hertzian model, equivalent K_n, K_s are
> (from PFC mannual):
>
> K_n = (2G sqrt(2R)/(3(1-nu)))sqrt(U_n);
>
> K_s = (2(G^2*3*(1-nu)R)^1/3/(2-nu))*F_n^1/3.
>
> So they are not constant but depend on the current overlap and sliding
> condition.
>
> Ning
>
> --
> 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
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

If it works for clumps, it's really chance.
I would guess it does not work, because for clumps the "l" branch should correspond to clump-clump distance not member-member (and the latest is probably what happens).
It's probbaly not very difficult to handle the clump case, anyone want to try?

Hi Bruno,
I need this feature (actually thanks for pointing this out) so I would be happy to code it. How do you suggest it?
Thanks, Chiara

On 23 Jun 2012, at 18:11, Chareyre wrote:

> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Chareyre proposed the following answer:
> If it works for clumps, it's really chance.
> I would guess it does not work, because for clumps the "l" branch should correspond to clump-clump distance not member-member (and the latest is probably what happens).
> It's probbaly not very difficult to handle the clump case, anyone want to try?
>
> --
> 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
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Something like this in getStress:

if (b1->clumpMember && b2->clumpMember){
l = O.bodies[b2->clumpId]->pos - O.bodies[b1->clumpId]->pos
} else {
//like now
}

It needs also to adapt this to clump vs. standalone.

For the periodic case, it needs to be shifted, but the current line should work without any change:
if (isPeriodic) branch-= scene->cell->hSize*I->cellDist.cast<Real>();

You can use triaxial test on clumps to validate (if boxes are used for boundaries, it gives an alternate calculation of stress). There is an example script from Janek IIRC.

On 23 Jun 2012, at 20:01, Chareyre wrote:

> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Chareyre proposed the following answer:
> Something like this in getStress:
>
> if (b1->clumpMember && b2->clumpMember){
> l = O.bodies[b2->clumpId]->pos - O.bodies[b1->clumpId]->pos
> } else {
> //like now
> }
OK, thank you.
>
> It needs also to adapt this to clump vs. standalone.
Yes, right.
>
> For the periodic case, it needs to be shifted, but the current line should work without any change:
> if (isPeriodic) branch-= scene->cell->hSize*I->cellDist.cast<Real>();
>
> You can use triaxial test on clumps to validate (if boxes are used for
> boundaries, it gives an alternate calculation of stress). There is an
> example script from Janek IIRC.
Yes, I know which script. So will try a comparison and let you know.
Chiara
>
> --
> 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
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

On 23 Jun 2012, at 20:01, Chareyre wrote:

> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Chareyre proposed the following answer:
> Something like this in getStress:
>
> if (b1->clumpMember && b2->clumpMember){
> l = O.bodies[b2->clumpId]->pos - O.bodies[b1->clumpId]->pos
> } else {
> //like now
> }
>
> It needs also to adapt this to clump vs. standalone.
>
> For the periodic case, it needs to be shifted, but the current line should work without any change:
> if (isPeriodic) branch-= scene->cell->hSize*I->cellDist.cast<Real>();
Bruno, why is it -=? In PeriIso*.cpp it used to be +=. Am I missing something?
Thank you, Chiara

>
> You can use triaxial test on clumps to validate (if boxes are used for
> boundaries, it gives an alternate calculation of stress). There is an
> example script from Janek IIRC.
>
> --
> 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
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Sure there is not another difference elsewhere giving the same result in
the end (e.g. sign of force)?

On 24 Jun 2012, at 11:15, Chareyre wrote:

> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Chareyre proposed the following answer:
> Sure there is not another difference elsewhere giving the same result in
> the end (e.g. sign of force)?
Yes, in fact there is no reversed sign inside the function getStress but instead there is in PeriIso*.cpp. The result should be the same in the end (Jan, have you done any test in these regards?). I gather it was decided not to keep reversing the sign of contact forces depending on the convention but I am not sure it is the best solution (I think it was a good reminder to not overlook that feature).

Chiara

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

>I think it was a good reminder to not overlook that feature

The "feature" was non-sense, so there is no need to remind it.
"-getStress()" makes more sense that "getStress(inverseSign=True)".

Else, we would need also sqrt(X,inverseSign=True), sin(X,inverseSign=True), etc. ;)

Hello Bruno,

I just now face a problem of the sign :-) why is branch defined as
id1 - id2
and not opposite?

Thanks
Jan

On 24.6.2012 19:01, Chareyre wrote:
> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Chareyre posted a new comment:
>> I think it was a good reminder to not overlook that feature
> The "feature" was non-sense, so there is no need to remind it.
> "-getStress()" makes more sense that "getStress(inverseSign=True)".
>
> Else, we would need also sqrt(X,inverseSign=True),
> sin(X,inverseSign=True), etc. ;)
>

I am doing some tests now, and the first result is that with minus

if (isPeriodic) branch -= scene->cell->hSize*I->cellDist.cast<Real>();

the results are not ok, while with

if (isPeriodic) branch += scene->cell->hSize*I->cellDist.cast<Real>();

the results are ok.. In the tests, I am using convention that positive
normal stress is tensile.. I will look on I->cellDist (I have never used
it) and let you know if I discover something :-)

My idea when thinking of not using reverseSign flag was that if you use
positive forces for tension, you would obtain positive stress for
tension and vice versa, but maybe this assumption is not correct..

Jan

On 24.6.2012 13:25, Chiara Modenese wrote:
> Question #199704 on Yade changed:
> https://answers.launchpad.net/yade/+question/199704
>
> Chiara Modenese proposed the following answer:
>
> On 24 Jun 2012, at 11:15, Chareyre wrote:
>
>> Question #199704 on Yade changed:
>> https://answers.launchpad.net/yade/+question/199704
>>
>> Chareyre proposed the following answer:
>> Sure there is not another difference elsewhere giving the same result in
>> the end (e.g. sign of force)?
> Yes, in fact there is no reversed sign inside the function getStress but instead there is in PeriIso*.cpp. The result should be the same in the end (Jan, have you done any test in these regards?). I gather it was decided not to keep reversing the sign of contact forces depending on the convention but I am not sure it is the best solution (I think it was a good reminder to not overlook that feature).
>
> Chiara
>
>
>> --
>> 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
>> Post to : <email address hidden>
>> Unsubscribe : https://launchpad.net/~yade-users
>> More help : https://help.launchpad.net/ListHelp

Jan,
Current branch makes tensile stress positive for ScGeom laws.
The problem in your case is (I guess) that you are using a Dem3Dof law, which invert the force convention (applied on b2 in ScGeom, applied on b1 in Dem3Dof). It's a long standing problem...

There was an agreement months ago to make ScGeom convention the default (hence, it needs to specify dem3Dof=True to make PeriTriaxController work correctly (https://github.com/yade/trunk/blob/master/pkg/dem/PeriIsoCompressor.hpp#L50)).

I would recommend to change the convention in Dem3Dof to avoid this, but it's up to you to decide since I think you are the only user/maintainer of associated code.
Currently you need -getStress() to get correct sign with Dem3Dof.

Hi,

I'm a little bit comfused following your discussion. It seems current version of YADE cannot calculate stress tensor and/or fabric tensor for clumps?

Ning

This question was expired because it remained in the 'Open' state without activity for the last 15 days.