yade-users team mailing list archive

[Nima Goudarzi <question652968@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Nima Goudarzi is still having a problem:
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

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.