yade-users team mailing list archive

[Anton Gladky <question253045@xxxxxxxxxxxxxxxxxxxxx>]

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
<question253045@xxxxxxxxxxxxxxxxxxxxx>:
> 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.