Question #253045 on Yade changed: https://answers.launchpad.net/yade/+question/253045
Anton Gladky posted a new comment: Hi Dominik, thanks again for the code analyze. I am glad, that the previous ViscElPM implementation gave correct results (except cs-error in the paper). I forgot, that kn1=2*kn... actually. I have committed changes [1], feel free to report, if you find something else. Could you please also prepare a verification script based on your example, which would compare analytical and numerical results explicitly? Like it is done here for example [2], but in a better way like in your case? I would add it to our check and test scripts, so the ViscEl model will be verified after each commit. [1] https://github.com/yade/trunk/commit/b1923146d8ffc0883f372b78e5ed1cba927dd801 [2] https://github.com/yade/trunk/blob/master/scripts/checks-and-tests/checks/checkViscElEng.py Thanks Anton 2014-08-26 16:02 GMT+02:00 Dominik Boemer <question253...@answers.launchpad.net>: > Question #253045 on Yade changed: > https://answers.launchpad.net/yade/+question/253045 > > Dominik Boemer posted a new comment: > Hello Anton, > > sorry for contacting you again. But I think that we missed something. > I have been reading [Pournin2001] and the source code again. Meanwhile > I have also been testing the code with the script you can find below. > What I found, is that the equivalent mass > > const Real massR = 2*mass1*mass2/(mass1+mass2); // WITH A FACTOR 2! > > delivers correct results. Why? Pournin's formulas are correct (except > the last expression in equation (22), as I mentionned earlier). This > means that there has to be another bug in the code. And this bug, is > the following: > > [This is currently in the source code.] > const Real massR = mass1*mass2/(mass1+mass2); > kn1 = kn2 = 1/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR; > cn1 = cn2 = -2.0 /Tc * log(En)*massR; > ks1 = ks2 = 2.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR; > cs1 = cs2 = -4.0/7.0 /Tc * log(Et)*massR; > > should be replaced by: > > [This should be in the source code.] > const Real massR = mass1*mass2/(mass1+mass2); > kn1 = kn2 = 2.0/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR; > cn1 = cn2 = -4.0 /Tc * log(En)*massR; > ks1 = ks2 = 4.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR; > cs1 = cs2 = -8.0/7.0 /Tc * log(Et)*massR; > > Why did I multiply kn1, kn2, cn1, cn2, ks1, ks2, cs1 and cs2 by 2 > instead of multiplying massR by 2 (as it was before)? Because massR is > given by the expression in [2001Pournin, just below equation (19)]. > kn1, kn2, cn1, ... is however not given explicitly by equation > (22)[2001Pournin]. This equation gives the resultant stiffness and > damping coefficients kn, ks, cn and cs (and not the individual > components kn1, kn2, ...). > > For instance, > > kn = 1.0/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR; > > according to [Pournin2001], and as kn1 = kn2, > > 1/kn = 1/kn1 + 1/kn2 <==> kn = kn1/2 <==> kn1 = 2*kn, which explains why > I multiplied the expressions of kn1, kn2, ... by 2. > > > To put it simple, we used the expressions of kn, ks, cn, ... in [Pournin2001] > for kn1, kn2, ks1, ... which is wrong, as kn1 = 2*kn, kn2 = 2*kn, ks1 = 2*ks > (when kn1 = kn2, ks1 = ks2, ...). > > As mentionned above, here is the script, which helped me to validate the > constitutive model in the normal direction. It shows that the contact > time is tc = 0.1 s and that the analytical velocity is -0.3 m/s when tc > = 0.1 s (this is the consequence of the normal restitution coefficient > en = 0.3). > > Best regards, > Dominik -- 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 : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp