← Back to team overview

yade-users team mailing list archive

Re: The critical time step?

 

On 4/5/08, yade-users-request@xxxxxxxxxxxxxxxx <
yade-users-request@xxxxxxxxxxxxxxxx> wrote:

> Send Yade-users mailing list submissions to
>        yade-users@xxxxxxxxxxxxxxxx
>
> To subscribe or unsubscribe via the World Wide Web, visit
>        https://lists.berlios.de/mailman/listinfo/yade-users
> or, via email, send a message with subject or body 'help' to
>        yade-users-request@xxxxxxxxxxxxxxxx
>
> You can reach the person managing the list at
>        yade-users-owner@xxxxxxxxxxxxxxxx
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of Yade-users digest..."
>
>
> Today's Topics:
>
>   1. The critical time step? (kan)
>   2. (no subject)
>   3. Re: The critical time step? (Bruno Chareyre)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Fri, 4 Apr 2008 10:40:41 -0500
> From: kan <nicessgg@xxxxxxxxx>
> Subject: [Yade-users] The critical time step?
> To: yade-users@xxxxxxxxxxxxxxxx
> Message-ID:
>        <45213f870804040840x267a75b4g3285992b589b1efb@xxxxxxxxxxxxxx>
> Content-Type: text/plain; charset="utf-8"
>
> Hey, guys,
>
> I have a problem to get the critical timestep, could you please give me
> hand?
> if I need to calculate the critical time step for each time steps, how can
> I
> calculate it?
>
> ------------------------------
>
> Message: 2
> Message-ID: <mailman.16.1207389622.4346.yade-users@xxxxxxxxxxxxxxxx>
>
> critical_time_step=2*sqrt(mass/k),
>
> my questions are:
> 1. is the mass the particle mass?
> 2. k here is kn or ks or some other value else?
> 3. I need to calculate all the critical_time_step for each particle, and
> then use the mininum one, right?
>
> and I hear that in the past, the critical_time_step need to calculated at
> each timestep, then why and how?
> if I do it at the beginning and get the mininum critical_time_step, and
> then
> keep it as the critical time step for the overall simulation, will it
> work?
>
> Thanks for your help.
>
> Kan.
> Message: 3
> Date: Fri, 04 Apr 2008 18:38:39 +0200
> From: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
> Subject: Re: [Yade-users] The critical time step?
> To: yade-users@xxxxxxxxxxxxxxxx
> Message-ID: <47F6598F.9030604@xxxxxxxxxxx>
> Content-Type: text/plain; charset=UTF-8; format=flowed
>
> kan a ?crit :
> > Hey, guys,
> >
> > I have a problem to get the critical timestep, could you please give
> > me hand?
> > if I need to calculate the critical time step for each time steps, how
> > can I calculate it?
> > From papers, I found that equation is:
> > critical_time_step=2*sqrt(mass/k),
> >
> > my questions are:
> > 1. is the mass the particle mass?
> Its the generalized mass. i.e. one mass for each degree of liberty :
> actual mass for 3 translation axis, and inertia matrix for the 3
> rotation axis. Its a bit difficult to uncouple rotational degrees of
> liberty because corresponding "mass" is a matrix...
> Cundall proposed to use diagonal terms as an approximation, its what I
> did in GlobalStiffnessTimestepper.


Thank you, Bruno.
I read the code, and found the code as following.

00055         Vector3r& stiffness = (static_cast<GlobalStiffness*>( ncb->
physicalActions<http://yade.berlios.de/doxygen/html/df/deb/classMetaBody.html#5b485c55dd71687385a781e046e0cf45>->find
(body->getId(),
globalStiffnessClassIndex<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#e54a0bd2675dce826ee3f06f636ce1f1>
).get()))->stiffness;
00056 Vector3r& Rstiffness = (static_cast<GlobalStiffness*>( ncb->
physicalActions<http://yade.berlios.de/doxygen/html/df/deb/classMetaBody.html#5b485c55dd71687385a781e046e0cf45>->find
(body->getId(),
globalStiffnessClassIndex<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#e54a0bd2675dce826ee3f06f636ce1f1>
).get()))->Rstiffness;

1. Then where is the stiffness and Rstiffness created, and how is it
created, I can not find the file. And what is the physical rules behind to
creat the global translational stiffness and global rotational stiffness?
thanks.

00083         Real dt = max( max (stiffness.X(), stiffness.Y()),
stiffness.Z() );
00084         if (dt!=0) {
00085                 dt =
sqrt(sdec->mass<http://yade.berlios.de/doxygen/html/d2/d85/classParticleParameters.html#c144beb4217eff38a9c0805c1b6ab5b2>/dt);
computedSomething<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#04f80a2ee4ce8d30f918178fe2a40af3>=
true;}

2. 1  here, it is sqrt(mass/stiffness), but I did not see like 2* as the
equation I have seen. why the 2 is lost or I am wrong?

00086         else dt = Mathr::MAX_REAL;
00087
00088         if (Rstiffness.X()!=0) {
00089                 dtx =
sdec->inertia<http://yade.berlios.de/doxygen/html/d7/d72/classRigidBodyParameters.html#8e165f286a97c3c1a7bfb0d4e8e703be>.X()/Rstiffness.X();
computedSomething<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#04f80a2ee4ce8d30f918178fe2a40af3>=
true;}
00090         else dtx = Mathr::MAX_REAL;
00091
00092         if (Rstiffness.Y()!=0) {
00093                 dty =
sdec->inertia<http://yade.berlios.de/doxygen/html/d7/d72/classRigidBodyParameters.html#8e165f286a97c3c1a7bfb0d4e8e703be>.Y()/Rstiffness.Y();
computedSomething<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#04f80a2ee4ce8d30f918178fe2a40af3>=
true;}
00094         else dty = Mathr::MAX_REAL;
00095
00096         if (Rstiffness.Z()!=0) {
00097                 dtz =
sdec->inertia<http://yade.berlios.de/doxygen/html/d7/d72/classRigidBodyParameters.html#8e165f286a97c3c1a7bfb0d4e8e703be>.Z()/Rstiffness.Z();
computedSomething<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#04f80a2ee4ce8d30f918178fe2a40af3>=
true;}
00098         else dtz = Mathr::MAX_REAL;
00099
00100
00101         //if (Rstiffness.X()!=0)
00102         //Real dtx = (Rstiffness.X()!=0 ?
sdec->inertia.X()/Rstiffness.X() : Mathr::MAX_REAL);
00103         //Real dty = (Rstiffness.Y()!=0 ?
sdec->inertia.Y()/Rstiffness.Y() : Mathr::MAX_REAL);
00104         //Real dtz = (Rstiffness.Z()!=0 ?
sdec->inertia.Z()/Rstiffness.Z() : Mathr::MAX_REAL);
00105         Real Rdt = sqrt( min( min (dtx, dty), dtz ) );

2. 2 the same question , why only sqrt( ), but not 2 or some other constant
factor.

00106
00107         //cerr << "Rstiffness.x()=" << Rstiffness.x() << "  " <<
Rstiffness.y() << "  " << Rstiffness.z() << endl;
00108         //cerr << "sdec->inertia=" << sdec->inertia.x() << " " <<
sdec->inertia.x() << " " << sdec->inertia.x() << endl;
00109         //cerr << "timesteps : dt=" << dt << " / Rdt=" << Rdt << endl;
00110
00111         dt =
0.709*timestepSafetyCoefficient<http://yade.berlios.de/doxygen/html/d6/d99/classGlobalStiffnessTimeStepper.html#63258ee4705252db25254e49900914dc>
*std::min(dt,Rdt);//0.709 = 1/sqrt(2)

3. here, what is the meaning of 0.709 where sqrt(2)/2, I understand the
"timestepSafetyCoefficient" is for safety purpose, but I cannot understand
the value of sqre(2)/2;




> > 2. k here is kn or ks or some other value else?
>
> Again generalized stiffness : a matrix including translational and
> rotational stiffnesses (taking into account all contacts on bodies). And
> again I used the diagonal terms to associate a scalar stiffness to each
> degree of liberty.


as question 1, how to calculate the matrix that includes translational and
rotational stiffnesses? or where is the function in Yade (or the file)?

 Thanks for your guildance very much.

Kan





> > 3. I need to calculate all the critical_time_step for each particle,
> > and then use the mininum one, right?
> yes
> >
> > and I hear that in the past, the critical_time_step need to calculated
> > at each timestep, then why and how?
> Contact created/deleted = stiffness changed.




> > if I do it at the beginning and get the mininum critical_time_step,
> > and then keep it as the critical time step for the overall simulation,
> > will it work?
> No, unless you set big security margin, which would mean a waste of time.
>
> Bruno
>
> --
>
> _______________
> Chareyre Bruno
> Maitre de conference
>
> Institut National Polytechnique de Grenoble
> Laboratoire 3S (Soils Solids Structures) - bureau E145
> BP 53 - 38041, Grenoble cedex 9 - France
> T?l : 33 4 56 52 86 21
> Fax : 33 4 76 82 70 43
> ________________
>
>
>
> ------------------------------
>
> _______________________________________________
> Yade-users mailing list
> Yade-users@xxxxxxxxxxxxxxxx
> https://lists.berlios.de/mailman/listinfo/yade-users
>
>
> End of Yade-users Digest, Vol 22, Issue 5
> *****************************************
>
_______________________________________________
Yade-users mailing list
Yade-users@xxxxxxxxxxxxxxxx
https://lists.berlios.de/mailman/listinfo/yade-users