Implementing Local Damping

Dear Nima, I think it is very sensible to avoid viscous damping for quasi-static problems, and I have good news for you:
1/ Cundall's damping is unrelated to the contact model, it appears only via the time integration scheme and it is already implemented.
2/ using the current HM with cn=cs=0 should not be a problem

All in all, you have nothing to implement except those little differences between current HM and the HM you have in mind.
You can do that directly by modifying HM law. If you can make it cleanly (backward compatible, ideally adding features without removing any, documented, examplified) we could consider including your changes in the main branch.
Regards
Bruno

Dear Bruno,

 Good to hear back from you and thanks so much for the reply. Actually I need some helps to make sure what I am implementing is true, Please forgive me for the lengthy email. I am sure it does not take you much time to take a look at my questions and correct me. The model, I'm implementing is not that far from the original HM model implemented by Chiara Modenese. I am trying to enrich rolling and twisting responses in this model (I will only use Law2_ScGeom_MindlinPhys_Mindlin as the basis of my own Law2). I have some basic questions, regarding what I am trying to implement. Although it seems too lengthy, my main concern is the possibility of passing some model parameters to Iphys (I think there is no need to define a new class of material but please correct me if I am wrong). The answers of questions 7,8 and 13 are criticl. The model (theoretically- I say theoretically because when the run the model, by considering local damping they set cn =cs=0 and subsequently cr = ct= 0 [these 2 are not free parameters]) considers 8 responses (4 contact responses and 4 damping responses). It has seven free parameters (P1 to P7) P1      E is the macro-scale modulus of elasticity. From which Kn is calculated as (Eq.1) Kn=2*r*E, where                                (Eq.2)r =(2ri*rj)/(ri+ri) is the common radius.(1)  ????: I have not found radius1 and radius2 as objects in any classes in YADE. Do I need to use lines 43 to 45 of HertzMindlin.cpp to implement (Eq.2):GenericSpheresContact* scg = YADE_CAST<GenericSpheresContact*>(interaction->geom.get());                      Real Da = scg->refR1>0 ? scg->refR1 : scg->refR2;                     Real Db = scg->refR2;

Or need to use radius1 and radius2 like ScGeom* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get() );                                               Real Da = geom->radius1;                                               Real Db = geom->radius2;
 In the original HM, the equivalent E is calculated from the young modulus of two interacting materials
(Real E = Ea*Eb/((1.-std::pow(Va,2))*Eb+(1.-std::pow(Vb,2))*Ea) (Line 54 of HertzMindlin.cpp)
(2)????: How my Macro-Scale E is implemented here? Do I still need to use  Real Ea=mat1->young;  and   (Line 35-HertzMindlin.cpp)                                                         Real Eb = mat2->young; (Line 36-HertzMindlin.cpp)
to construct a macro-scale E or there is another way to do so (For example defining a unique E for the whole material, if this is possible). Currently, I have implemented the normal stiffness as    Kno=4*Ea*Da*Eb*Db/(Ea*Da+Eb*Db) (line 57 of  original Hertz Mindlin.cpp will be this)   which only will be equal to (1) if Ea=Eb In this model (3)Ks=Kn/ξ where  P2      ξ is the ratio of the normal to the tangential contact stiffness. I have declared ξ in Ip2 as

       ((Real,Xi,1.0,, "Ratio of the normal to the tangential contact stiffness"))

 (3)???? Is this a correct place for declaring ξ or I need to define a new class (for example a new material class) to declare this constant?
If I declare ξ in Ip2 Kso=Kno/Xi (line 58 of  original Hertz Mindlin.cpp will be this)   P3     β is a shape parameter used to consider the effects of particle shape on the overall mechanical behavior of granular materials. It, actually,  links the contact radius Ṝ to common radius (Eq.2) as

 (Eq.3) Ṝ=βr

I, again, have declared β in Ip2 as

       ((Real,beta,0.0,, "Dimensionless shape parameter linking the contact radius R bar to the r"))

 (4)???? Is this a correct place for declaring β?(Actually, I think, there is no need to write a new class of material).    This  Ṝ is so important to me. It will be used in many places throughout the code in the calculation of rotational and twisting  stiffnesses and rolling and twisting damping coefficients (see 6 and 7 below): In the original HM, Kr  and Kt are input parameters (I think they should be given by user)   In the model, I am implementing, Kr  and Kt are not free parameters   (Eq.4) Kr=0.25KnṜ2

 (Eq.5) Kt=0.5KsṜ2

 If I declare both  ξ and β in Ip2   Kr=0.25*Kno*std::pow((beta*r),2); (line 70 of  original Hertz Mindlin.cpp will be this)
 (5) (a c++ question) ???? Do I need to use std::pow or a simple pow does the job?
Kt=0.5*(Kno/Xi)*std::pow((beta*r),2); (line 71 of  original Hertz Mindlin.cpp will be this)
(6) ???? Are these equations true (Regarding that I have declared both ξ and β in Ip2?
(7) ???? Can I calculate Ṝ and pass it to my Iphys using contactPhysics pointer?  I later will need to use this Ṝ for calculation of rolling and twisting damping coefficients. contactPhysics->R bar =r*beta
Real r =2*((geom->radius1*geom->radius2)/(geom->radius1*geom->radius2)
or                                        Real r =2*((scg->refR1*scg->refR2)/(geom->refR1*geom->refR2) (If I follow original HM code)

I am doubtful of passing R bar to Iphys, because I know that Iphys is the home of non-geometrical contact properties but Ṝ is somehow related to the geometry of contact. If no what do I need to do to be able to call Ṝ (or even beta*r) later in the calculation of rolling and twisting damping coefficients and peak resistances in rolling and twisting directions?
P4       ξc is a local crushing parameter describing the effects of local asperity crushing and related to the hardness of the particle mineral material. It is introduced to control the peak rolling resistance. (See c below)Actually, the model imposes a Mohr-Columb like criterion in rolling and twisting directions to introduce the peak resistances in each direction (Like Mohr-Columb criterion in tangential direction which has already been implemented in original HM. See b above)

(8) ???? In which class do I need to implement ξc? Do I need a new class of material for this? or I can declare it in Ip2 and pass it to Iphys ?
In Ip2 (.hpp)
      ((Real,Xi c,2.1,, "Local crushing parameter describing the effects of local asperity crushing"))
In Iphys(.hpp)      ((Real,XI c,2.1,, "Local crushing parameter describing the effects of local asperity crushing"))
In .cpp file
contactPhysics->XI c =Xi c s;
I will later call XI c using a pointer (phys for example) to implement a Mohr-Columb like rolling peak resistance
 Real maxMr =0.25*phys-> XI c*sdec->R bar* Fn(sdec is a pointer to R bar in Iphys) This is only true if I can pass r*beta to R bar and Xi c to XI c in my Iphys. (subjects of questions 7 and 8)
and
 Real maxMt =0.25*phys-> tangensOfFrictionAngle*sdec->R bar* Fn
(9) ???? I plan to implement the codes for peak resistances immediately after  /* MOHR-COULOMB law */(lines 389 to 421 of Original HM.cpp). As the peak resistances in rolling and twisting directions require to use Fn (see c&d above), do I need to update Fn  again regarding this fact that Fn , has been updated once in/* VISCOUS DAMPING (Normal direction) */ (lines 337 to 366 of Original HM.cpp)
(10)???? If I implement the peak resistances in rolling and twisting directions, I, think, I don't need to implement anything regarding plasticity condition in either direction. Do you think I am right?
P5         µ is the inter-particle friction coefficient (This has already been implemented)
P6&P7         cn and cs  are normal and shear damping coefficients.

In the model, I am implementing, these two not only are used in the calculation of normal and tangential damping forces but also are used for calculation of cr and ct these(rolling and twisting damping coefficients) as :(Eq.5) cr=0.25cnṜ2

(Eq.6) ct=0.5csṜ2
a) If the simulations are quasi-static, local damping has some advantages over viscous contact damping. This is what the paper has proposed
Note that viscous contact damping and local damping (PFC3Dmanual [53] has details on this form of damping) have no significanteffects in quasi-static behavior. Since local damping have severaladvantages identified in PFC3D [53], this study therefore used thelocal damping. A critical damping coefficient of 0.7 was chosen bytrial and error so that the particulate system is properly-damped,i.e., neither over-damped nor under-damped. The viscous contactdamping was turned off as cn = cs = 0. Viscous contact dampingcan be the proper damping in future simulations of dynamic problems.
This was a question I asked you yesterday and you answered that using the current HM with cn = cs=0 should not be a problem. If you refer to the original HM (/* DAMPING COEFFICIENTS */Lines 297 to 317), Damping coefficients have been defined for both linear and nonlinear elastic cases

(11)???? Do I need damping for both linear and nonlinear cases?

(12)???? If the answer of question 11 is yes, do I need to erase whatever related to Cn-crit and Cs-crit in linear damping and the formula (which calculates cn using alpha) in nonlinear damping since I need to directly set cs = cs=0?
In this case (all lines related to critical damping will be removed from the code)
                                                                                   Real cn=0;
                                                                                   Real cs=0;
Is this true?
(13) If I don't need to implement anything for numerical (local) damping in my code, How A critical damping coefficient of 0.7 (in the paper above) can be assigned (this is the default local damping value in PFC3D prior to ver.5). Will this be a part of writing the script?

b) I prefer to write a code which has the flexibility to account for dynamic simulations also. In this case I need to include viscous damping coefficients in all 4 directions:
(14)???? Do I need dampings for both linear and nonlinear cases in this situation?

(15)???? If the answer of question 14 is yes, can I immediately implement the equations for cr and ct after calculating cn and cs?

Something like this (Please take a look at lines in red)
|
|
|         /************************/ |
|         /* DAMPING COEFFICIENTS */ |
|         /************************/ |
|         |
|         // Inclusion of local damping if requested |
|         // viscous damping is defined for both linear and non-linear elastic case |
|         if (useDamping && LinDamp){ |
|                 Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1-                                         >isDynamic()) ? de1->mass : (de1->                 mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set                                 it equal to the one of the                  dynamic body |
|                 //Real mbar = de1->mass*de2->mass / (de1->mass + de2->mass); // equivalent mass |
|                 Real Cn_crit = 2.*sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction) |
|                 Real Cs_crit = 2.*sqrt(mbar*phys->ks); // Critical damping coefficient (shear direction) |
|                 // Note: to compare with the analytical solution you provide cn and cs directly (since here we used a different                   method to define c_crit) |
|                 cn = Cn_crit*phys->betan; // Damping normal coefficient |
|                 cs = Cs_crit*phys->betas; // Damping tangential coefficient
|                cr = 0.25*phys->R bar*cn // Damping rolling coefficient               ct = 0.5*phys->R bar*cs // Damping twisting coefficient (True???)
 |

 |
|                 if(phys->kn<0 || phys->ks<0){ cerr<<"Negative stiffness kn="<<phys->kn<<" ks="<<phys->ks<<" for ##"<<b1->getId()                  <<"+"<<b2->getId()<<", step "<<scene->iter<<endl; } |
|         } |
|         else if (useDamping){ // (see Tsuji, 1992) |
|                 Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ?                    de1->mass : (de1->                  mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to                 the one of the                  dynamic body |
|                 cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient, see also [Antypov2011] eq. 10
                For me (16)???? I  need to call coefficients of restitution directly (not calculating alpha and passing it to Iphys). Can I pass en and es from Ip2 to Iphys and use it here?
In Ip2 (.hpp) [This already exists in original HM]
|                         ((shared_ptr<MatchMaker>,en,,,"Normal coefficient of restitution $e_n$.")) |
|                         ((shared_ptr<MatchMaker>,es,,,"Shear coefficient of restitution $e_s$.")) |

In Iphys(.hpp)
(17)???? If passing from Ip2 to Iphys is possible, How I can declare en and es in Iphys?

In .cpp file (here) instead of  cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient, see also [Antypov2011] eq. 10
cn = (2*sqrt(mbar*kn)*ln phys->en [this is planned to be addeed in Iphys])/(sqrt(pow(ln phys->en,2)+pow(Mathr::PI,2)))  // normal viscous coefficient |
|  cs = cn; // same value for shear viscous coefficient(18)???? Is this true? Do I need to use a separate formula for cs?cs = (I will write my own cs)
cr = 0.25*phys->R bar*cn // Damping rolling coefficientct = 0.5*phys->R bar*cs // Damping twisting coefficient (True???)
 |
|         } |

 If viscous damping will be used (as in dynamic simulations), I will activate damping only in elastic regimes for tangential, rolling and twisting directions (as implemented for tangential direction in original HM). Following is the summary for derived formulations I am trying to implement |

 |

I have some other minor issues in the code, I'll ask after getting responses for asked questions.

Sincerely Yours,
Nima

>my main concern is the possibility of passing some model parameters to Iphys (I think there is no need to define a new class of material but please correct me if I am wrong)

It is possible and indeed simpler usually. [*] Adding the parameters to a material class gives more flexibility when mixing particles of different properties.

(1) radius1 is an alias of refR1, you can use both names equivalently (ScGeom.hpp:58)

(2) I see no need to change current code since it does what you need as a special case.

(3) ξ can be where you have it now or in the material class, see [*]

(4) same answer. In this case it is obvious that if you want to have particles of different shapes (different β) you need to have it in the material class, else all shapes are the same as defined by the Ip2.

(5) if both versions compile they probably do the same thing

(6) I don't see a problem

(7) yes (is it different from the assignement of e.g. contactPhysics->kno?)

(8) If ξc is like a damage parameter with different evolutions at each contact then it needs to be IPhys::ξc indeed. You don't "declare" it in IP2, you may only assign it there. But even the assignement seems to be in Iphys(.hpp) a.t.m.:
((Real,XI c,2.1,, "Local crushing parameter ..."))

(9) This question sounds like "which is the right constitutive model?". I don't know. This will be your modeling assumption. :)
If you plan to keep viscous damping at 0 it should not be a direct problem for you though.

(10) I would say yes, but I'm not sure about the question. "Plasticity condition" just means to keep forces/torques bounded by some max values.

(11) I don't know.

(12) No need to erase anything since those lines are conditional already (l.303,313), your changes should correspond to an additional option.

(13) Newton::damping is the value of Cundall's/PFC "non-viscous" damping (the one with default 0.7 in PFC). Again, it has nothing to do with the contact model (it is thus irrelevant to distinguish linear vs. non-linear models). Newton::damping is set in the scripts.

(14) depends on you

(15) why not? :)

Wow, I made it. ;)
Bruno

Dear Bruno,

I highly appreciate your kind attention toward my request. Many of your answers are extremely helpful. I am currently working on writing the code but still have some issues in energy terms and displacements (I am trying to adapt what have been implemented for shear to rolling and twisting but have some difficulties). I will be back soon to ask some extra questions regarding my new issues.

Sincerely yours,

Nima

> On Jul 26, 2017, at 7:53 AM, Bruno Chareyre <email address hidden> wrote:
>
> Your question #652968 on Yade changed:
> https://answers.launchpad.net/yade/+question/652968
>
> Status: Open => Answered
>
> Bruno Chareyre proposed the following answer:
>> my main concern is the possibility of passing some model parameters to
> Iphys (I think there is no need to define a new class of material but
> please correct me if I am wrong)
>
> It is possible and indeed simpler usually. [*] Adding the parameters to
> a material class gives more flexibility when mixing particles of
> different properties.
>
> (1) radius1 is an alias of refR1, you can use both names equivalently
> (ScGeom.hpp:58)
>
> (2) I see no need to change current code since it does what you need as
> a special case.
>
> (3) ξ can be where you have it now or in the material class, see [*]
>
> (4) same answer. In this case it is obvious that if you want to have
> particles of different shapes (different β) you need to have it in the
> material class, else all shapes are the same as defined by the Ip2.
>
> (5) if both versions compile they probably do the same thing
>
> (6) I don't see a problem
>
> (7) yes (is it different from the assignement of e.g.
> contactPhysics->kno?)
>
> (8) If ξc is like a damage parameter with different evolutions at each contact then it needs to be IPhys::ξc indeed. You don't "declare" it in IP2, you may only assign it there. But even the assignement seems to be in Iphys(.hpp) a.t.m.:
> ((Real,XI c,2.1,, "Local crushing parameter ..."))
>
> (9) This question sounds like "which is the right constitutive model?". I don't know. This will be your modeling assumption. :)
> If you plan to keep viscous damping at 0 it should not be a direct problem for you though.
>
> (10) I would say yes, but I'm not sure about the question. "Plasticity
> condition" just means to keep forces/torques bounded by some max values.
>
> (11) I don't know.
>
> (12) No need to erase anything since those lines are conditional already
> (l.303,313), your changes should correspond to an additional option.
>
> (13) Newton::damping is the value of Cundall's/PFC "non-viscous" damping
> (the one with default 0.7 in PFC). Again, it has nothing to do with the
> contact model (it is thus irrelevant to distinguish linear vs. non-
> linear models). Newton::damping is set in the scripts.
>
> (14) depends on you
>
> (15) why not? :)
>
> Wow, I made it. ;)
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/652968/+confirm?answer_id=2
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/652968
>
> You received this question notification because you asked the question.

Dear Bruno,
Thanks so much again and sorry for rising up this again. I am still a little confused with applying local damping:
a very quick question on your comment
(13) Newton::damping is the value of Cundall's/PFC "non-viscous" damping
> (the one with default 0.7 in PFC). Again, it has nothing to do with the
> contact model (it is thus irrelevant to distinguish linear vs. non-
> linear models). Newton::damping is set in the scripts.
                 /************************/

               /* DAMPING COEFFICIENTS */                    FROM HM.cpp

               /************************/

               // Inclusion of local damping if requested

               // viscous damping is defined for both linear and non-linear elastic case

               if (useDamping && LinDamp){

                              Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic()

                                 && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2-

                                 >mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of

                                 the dynamic body

                              //Real mbar = de1->mass*de2->mass / (de1->mass + de2->mass); // equivalent mass

                              Real Cn_crit = 2.*sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction)

                              Real Cs_crit = 2.*sqrt(mbar*phys->ks); // Critical damping coefficient (shear direction)

                              // Note: to compare with the analytical solution you provide cn and cs directly (since here

                                we used a different method to define c_crit)

                              cn = Cn_crit*phys->betan; // Damping normal coefficient

                              cs = Cs_crit*phys->betas; // Damping tangential coefficient

                              if(phys->kn<0 || phys->ks<0){ cerr<<"Negative stiffness kn="<<phys->kn<<"

                                ks="<<phys->ks<<" for ##"<<b1->getId()<<"+"<<b2->getId()<<", step "<<scene-

                                >iter<<endl; }

               }

               else if (useDamping){ // (see Tsuji, 1992)

                              Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic()

                                 && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2-

                                  >mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of

                                 the dynamic body

                              cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient, see also

                                [Antypov2011] eq. 10

                              cs = cn; // same value for shear viscous coefficient

               }

Do you mean I should keep unchanged everything included in the above code and set the local damping in the script? or I need to set cn=cs=0 wherever I find them in the code? If possible, would you please apply the required changes (if needed) in the above code if possible?

Sincerely yours,
Nima

    On Thursday, July 27, 2017 4:57 PM, Nima Goudarzi <email address hidden> wrote:

 Your question #652968 on Yade changed:
https://answers.launchpad.net/yade/+question/652968

    Status: Answered => Open

You are still having a problem:
Dear Bruno,

I highly appreciate your kind attention toward my request. Many of your
answers are extremely helpful. I am currently working on writing the
code but still have some issues in energy terms and displacements (I am
trying to adapt what have been implemented for shear to rolling and
twisting but have some difficulties). I will be back soon to ask some
extra questions regarding my new issues.

Sincerely yours,

Nima

> On Jul 26, 2017, at 7:53 AM, Bruno Chareyre <email address hidden> wrote:
>
> Your question #652968 on Yade changed:
> https://answers.launchpad.net/yade/+question/652968
>
>    Status: Open => Answered
>
> Bruno Chareyre proposed the following answer:
>> my main concern is the possibility of passing some model parameters to
> Iphys (I think there is no need to define a new class of material but
> please correct me if I am wrong)
>
> It is possible and indeed simpler usually. [*] Adding the parameters to
> a material class gives more flexibility when mixing particles of
> different properties.
>
> (1) radius1 is an alias of refR1, you can use both names equivalently
> (ScGeom.hpp:58)
>
> (2) I see no need to change current code since it does what you need as
> a special case.
>
> (3) ξ can be where you have it now or in the material class, see [*]
>
> (4) same answer. In this case it is obvious that if you want to have
> particles of different shapes (different β) you need to have it in the
> material class, else all shapes are the same as defined by the Ip2.
>
> (5) if both versions compile they probably do the same thing
>
> (6) I don't see a problem
>
> (7) yes (is it different from the assignement of e.g.
> contactPhysics->kno?)
>
> (8) If ξc is like a damage parameter with different evolutions at each contact then it needs to be IPhys::ξc indeed. You don't "declare" it in IP2, you may only assign it there. But even the assignement seems to be in Iphys(.hpp) a.t.m.:
> ((Real,XI c,2.1,, "Local crushing parameter ..."))
>
> (9) This question sounds like "which is the right constitutive model?". I don't know. This will be your modeling assumption. :)
> If you plan to keep viscous damping at 0 it should not be a direct problem for you though.
>
> (10) I would say yes, but I'm not sure about the question. "Plasticity
> condition" just means to keep forces/torques bounded by some max values.
>
> (11) I don't know.
>
> (12) No need to erase anything since those lines are conditional already
> (l.303,313), your changes should correspond to an additional option.
>
> (13) Newton::damping is the value of Cundall's/PFC "non-viscous" damping
> (the one with default 0.7 in PFC). Again, it has nothing to do with the
> contact model (it is thus irrelevant to distinguish linear vs. non-
> linear models). Newton::damping is set in the scripts.
>
> (14) depends on you
>
> (15) why not? :)
>
> Wow, I made it. ;)
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/652968/+confirm?answer_id=2
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/652968
>
> You received this question notification because you asked the question.

--
You received this question notification because you asked the question.

Hi,
This code you mention is related to viscous contact damping, not Cundall's damping. It appears to be doing nothing if "useDamping" is turned false, so I guess you can do just that and then there is no need to worry about cn, cs.
I mean to set useDamping=False in simulation scripts, there is absolutely nothing to change in the source code as I see it.
Cheers
Bruno

 Dear Bruno,
Thanks so much for all your helps. I finally implemented my model in YADE but got some new problems in modelling uniaxial compaction in the python script. I am trying to compact a soil sample layer by layer to reach a target void ratio for each layer and finally a target void ratio for the whole assembly. Worthy to say that, I'll have two distinct scripts. One for uniaxial compaction and one for triaxial test that is to say that I'll export final positions of spheres (after compaction) to my triaxial script. I am familiar with triaxStressController but prefer not to use it since I don't want to use the lubrication concept and expanding the particle size (or shrinking the walls). My final test is a trixal test but I prefer to use a very simple method for compacting using a plate from top of each layer. I mean, I want to use a facetBox as a mold and a simple plate. I have modeled the gravity deposition of the first layer and also put the plate on top of the layer and applied some constant velocity in z direction for moving the plate. Everything goes well but I have a problem with stopping the plate when the porosity of the layer reaches my target porosity. I got some hints from some YADE users regarding this to for checking the voxelPorosity  against my target porosity to stop the plate. Here are my questions:
1- I don't know how I can check the updating porosity against my target porosity since the volume of my layer is changing due to the compaction and as far as I know, voxelPorosity needs a height vector (let's say mx) which in my case is not constant and is gradually deceasing. Is there a way to stop the plate with target porosity in my case or I am enforced to use triaxStressController?
2- My parameters are updating from compaction to triaxial (due to numerical reasons). One of the reasons I prefer to use two scripts is this. Actually, I don't know if there is an approach to update material parameters in one single script. Please give me some hints if there is a solution.

Here is a draft of my script for compacting the first layer. It doesn't work from yade.pack import *from yade import utilsfrom yade import exportfrom yade import plotn_s = 2000n_band = 26# For now I'm using Sobieski's PSD. I need to change it to my own PSD latertargetPorosity= 0.422##################################################   Particle Size Distribution   ##################################################psdSizes=[5.8598400e+00,5.8761600e+00,5.8924800e+00,5.9088000e+00,5.9251200e+00,5.9414400e+00,5.9577600e+00,5.9740800e+00,5.9904000e+00,6.0067200e+00,6.0230400e+00,6.0393600e+00,6.0556800e+00,6.0720000e+00,6.0883200e+00,6.1046400e+00,6.1209600e+00,6.1372800e+00,6.1536000e+00,6.1699200e+00,6.1862400e+00,6.2025600e+00,6.2188800e+00,6.2352000e+00,6.2515200e+00,6.2678400e+00]psdCumm=[0.0,8.4945735e-05,3.5804112e-04,1.1512195e-03,3.2324246e-03,8.1658647e-03,1.8731095e-02,3.9172232e-02,7.4902028e-02,1.3132521e-01,2.1182373e-01,3.1558203e-01,4.3640887e-01,5.6352779e-01,6.8435463e-01,7.8811293e-01,8.6861145e-01,9.2503463e-01,9.6076443e-01,9.8120556e-01,9.9177079e-01,9.9670423e-01,9.9878544e-01,9.9957862e-01,9.9985171e-01,1.0]
i = 0while i < n_band:    psdSizes[i]=psdSizes[i]/1000. #scaling from [mm] to [m]    print "\r psdSizes:",i,psdSizes[i]    i = i +1
#If the bed volume is defined, then the average diameter of particles (when the target porosity is reached) is the same like in the measurement###########################################################   facetBox   ############################################################ Compaction will be in Y-Directionl = 3.5e-3h = 7.1e-3vol = l*h*l# You maybe don't need this if you use compaction from the top of the sample (using plate)
mn = Vector3(0,0,0)mx = Vector3(l,l,h/5)
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))O.bodies.append(geom.facetBox((l/2,l/2,h/2),(l/2,l/2,h/2),wallMask=31,material='walls'))# This will be an open-top boxsp = yade._packSpheres.SpherePack()sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=n_s,distributeMass=1)O.bodies.append([sphere(s[0],s[1],material='spheres') for s in sp])Porosity = voxelPorosity(200,mn,mx)print "\r Porosity: ",Porosity,O.engines=[   ForceResetter(),   # sphere, facet, wall   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),   InteractionLoop(      # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],      [Ip2_FrictMat_FrictMat_FrictPhys()],      [Law2_ScGeom_FrictPhys_CundallStrack()]   ),   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),   # the label creates an automatic variable referring to this engine   # we use it below to change its attributes from the functions called   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),]O.dt=.5*PWaveTimeStep()
def checkUnbalanced():   # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable   if O.iter<5000: return    # the rest will be run only if unbalanced is < .1 (stabilized packing)   if unbalancedForce()>.1: return    # add plate at the position on the top of the packing   # the maximum finds the z-coordinate of the top of the topmost particle   O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))   global plate        # without this line, the plate variable would only exist inside this function   plate=O.bodies[-1]  # the last particles is the plate   # Wall objects are "fixed" by default, i.e. not subject to forces   # prescribing a velocity will therefore make it move at constant velocity (downwards)   plate.state.vel=(0,0,-.1)   # start plotting the data now, it was not interesting before   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]   # next time, do not call this function anymore, but the next one (unloadPlate) instead   checker.command='unloadPlate()'def unloadPlate():   # if the porosity of the layer goes below the target porosity stop unloading   maxz=max([b.state.pos[2]+b.shape.radius for b in O.bodies if is instance(b.shape,Sphere)])   mx=maxz# Iknow this doesn't work   Por_osity=voxelPorosity(200,mn,mx)   if Por_osity <targetPorosity:      plate.state.vel=(0,0,0)

Thanks so much,
Sincerely yours
Nima On Monday, July 31, 2017, 7:27:43 AM CDT, Bruno Chareyre <email address hidden> wrote:

 Your question #652968 on Yade changed:
https://answers.launchpad.net/yade/+question/652968

    Status: Open => Answered

Bruno Chareyre proposed the following answer:
Hi,
This code you mention is related to viscous contact damping, not Cundall's damping. It appears to be doing nothing if "useDamping" is turned false, so I guess you can do just that and then there is no need to worry about cn, cs.
I mean to set useDamping=False in simulation scripts, there is absolutely nothing to change in the source code as I see it.
Cheers
Bruno

--
If this answers your question, please go to the following page to let us
know that it is solved:
https://answers.launchpad.net/yade/+question/652968/+confirm?answer_id=5

If you still need help, you can reply to this email or go to the
following page to enter your feedback:
https://answers.launchpad.net/yade/+question/652968

You received this question notification because you asked the question.

 Excuse me Bruno. Attached please find the script. The last one (in the email) is not readable I think.
Sincerely yours,
Nima On Tuesday, October 17, 2017, 5:10:16 PM CDT, Nima Goudarzi <email address hidden> wrote:

 Your question #652968 on Yade changed:
https://answers.launchpad.net/yade/+question/652968

    Status: Answered => Open

You are still having a problem:
 Dear Bruno,
Thanks so much for all your helps. I finally implemented my model in YADE but got some new problems in modelling uniaxial compaction in the python script. I am trying to compact a soil sample layer by layer to reach a target void ratio for each layer and finally a target void ratio for the whole assembly. Worthy to say that, I'll have two distinct scripts. One for uniaxial compaction and one for triaxial test that is to say that I'll export final positions of spheres (after compaction) to my triaxial script. I am familiar with triaxStressController but prefer not to use it since I don't want to use the lubrication concept and expanding the particle size (or shrinking the walls). My final test is a trixal test but I prefer to use a very simple method for compacting using a plate from top of each layer. I mean, I want to use a facetBox as a mold and a simple plate. I have modeled the gravity deposition of the first layer and also put the plate on top of the layer and applied some constant velocity in z direction for moving the plate. Everything goes well but I have a problem with stopping the plate when the porosity of the layer reaches my target porosity. I got some hints from some YADE users regarding this to for checking the voxelPorosity  against my target porosity to stop the plate. Here are my questions:
1- I don't know how I can check the updating porosity against my target porosity since the volume of my layer is changing due to the compaction and as far as I know, voxelPorosity needs a height vector (let's say mx) which in my case is not constant and is gradually deceasing. Is there a way to stop the plate with target porosity in my case or I am enforced to use triaxStressController?
2- My parameters are updating from compaction to triaxial (due to numerical reasons). One of the reasons I prefer to use two scripts is this. Actually, I don't know if there is an approach to update material parameters in one single script. Please give me some hints if there is a solution.

Here is a draft of my script for compacting the first layer. It doesn't work from yade.pack import *from yade import utilsfrom yade import exportfrom yade import plotn_s = 2000n_band = 26# For now I'm using Sobieski's PSD. I need to change it to my own PSD latertargetPorosity= 0.422##################################################   Particle Size Distribution   ##################################################psdSizes=[5.8598400e+00,5.8761600e+00,5.8924800e+00,5.9088000e+00,5.9251200e+00,5.9414400e+00,5.9577600e+00,5.9740800e+00,5.9904000e+00,6.0067200e+00,6.0230400e+00,6.0393600e+00,6.0556800e+00,6.0720000e+00,6.0883200e+00,6.1046400e+00,6.1209600e+00,6.1372800e+00,6.1536000e+00,6.1699200e+00,6.1862400e+00,6.2025600e+00,6.2188800e+00,6.2352000e+00,6.2515200e+00,6.2678400e+00]psdCumm=[0.0,8.4945735e-05,3.5804112e-04,1.1512195e-03,3.2324246e-03,8.1658647e-03,1.8731095e-02,3.9172232e-02,7.4902028e-02,1.3132521e-01,2.1182373e-01,3.1558203e-01,4.3640887e-01,5.6352779e-01,6.8435463e-01,7.8811293e-01,8.6861145e-01,9.2503463e-01,9.6076443e-01,9.8120556e-01,9.9177079e-01,9.9670423e-01,9.9878544e-01,9.9957862e-01,9.9985171e-01,1.0]
i = 0while i < n_band:    psdSizes[i]=psdSizes[i]/1000. #scaling from [mm] to [m]    print "\r psdSizes:",i,psdSizes[i]    i = i +1
#If the bed volume is defined, then the average diameter of particles (when the target porosity is reached) is the same like in the measurement###########################################################   facetBox   ############################################################ Compaction will be in Y-Directionl = 3.5e-3h = 7.1e-3vol = l*h*l# You maybe don't need this if you use compaction from the top of the sample (using plate)
mn = Vector3(0,0,0)mx = Vector3(l,l,h/5)
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))O.bodies.append(geom.facetBox((l/2,l/2,h/2),(l/2,l/2,h/2),wallMask=31,material='walls'))# This will be an open-top boxsp = yade._packSpheres.SpherePack()sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=n_s,distributeMass=1)O.bodies.append([sphere(s[0],s[1],material='spheres') for s in sp])Porosity = voxelPorosity(200,mn,mx)print "\r Porosity: ",Porosity,O.engines=[   ForceResetter(),   # sphere, facet, wall   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),   InteractionLoop(      # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],      [Ip2_FrictMat_FrictMat_FrictPhys()],      [Law2_ScGeom_FrictPhys_CundallStrack()]   ),   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),   # the label creates an automatic variable referring to this engine   # we use it below to change its attributes from the functions called   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),]O.dt=.5*PWaveTimeStep()
def checkUnbalanced():   # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable   if O.iter<5000: return    # the rest will be run only if unbalanced is < .1 (stabilized packing)   if unbalancedForce()>.1: return    # add plate at the position on the top of the packing   # the maximum finds the z-coordinate of the top of the topmost particle   O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))   global plate        # without this line, the plate variable would only exist inside this function   plate=O.bodies[-1]  # the last particles is the plate   # Wall objects are "fixed" by default, i.e. not subject to forces   # prescribing a velocity will therefore make it move at constant velocity (downwards)   plate.state.vel=(0,0,-.1)   # start plotting the data now, it was not interesting before   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]   # next time, do not call this function anymore, but the next one (unloadPlate) instead   checker.command='unloadPlate()'def unloadPlate():   # if the porosity of the layer goes below the target porosity stop unloading   maxz=max([b.state.pos[2]+b.shape.radius for b in O.bodies if is instance(b.shape,Sphere)])   mx=maxz# Iknow this doesn't work   Por_osity=voxelPorosity(200,mn,mx)   if Por_osity <targetPorosity:      plate.state.vel=(0,0,0)

Thanks so much,
Sincerely yours
Nima    On Monday, July 31, 2017, 7:27:43 AM CDT, Bruno Chareyre <email address hidden> wrote:

 Your question #652968 on Yade changed:
https://answers.launchpad.net/yade/+question/652968

    Status: Open => Answered

Bruno Chareyre proposed the following answer:
Hi,
This code you mention is related to viscous contact damping, not Cundall's damping. It appears to be doing nothing if "useDamping" is turned false, so I guess you can do just that and then there is no need to worry about cn, cs.
I mean to set useDamping=False in simulation scripts, there is absolutely nothing to change in the source code as I see it.
Cheers
Bruno

--
If this answers your question, please go to the following page to let us
know that it is solved:
https://answers.launchpad.net/yade/+question/652968/+confirm?answer_id=5

If you still need help, you can reply to this email or go to the
following page to enter your feedback:
https://answers.launchpad.net/yade/+question/652968

You received this question notification because you asked the question.

You received this question notification because you asked the question.

1- Is your question "how can I calculate porosity?". I would think it is trivial if the volume is defined by facets (you can use their positions).
2- There is no problem in changing parameters of O.materials. You'll need to make sure existing interactions are updated though (possibly looping over them to update some values).

Your copy-paste failed somehow, so I'm unable to read your script without line-breaks. I would suggest pasting it in the launchpad interface (I guess you sent by email?).

I must say I would be more inclined to help someone trying to avoid export/import of positions, because it will cause a number of superfluous problems later.
Positions alone do not define a mechanical state. Most of the information is lost in the pipe.

Bruno