equivalent damping calculation in ViscElPhys

>if one damping parameter of the contact is zero, then as the two dashpot are considered in series, the equivalent damping should be zero.

I tend to fully agree if you speak of single dashpots in series. However, if you say each particle's surface is deformable in a visco-elastic manner and the deformations of the two sufaces are cumulative (reason behind the harmonic average of stiffness) then it could make more sense to say that you have two parallel spring-dashpots systems in series, like model B in [1] except that Mt=0.

Then the question becomes: what is the equivalent viscosity (if possible to define one, it may not be the case[*]) of such system? I don't think it is cn = cn1cn2/(cn1+cn2). It should not be zero either when one of the particles is purely elastic (cn=0), since the second one is still viscous.

[1] http://www.sccs.swarthmore.edu/users/06/dluong1/e12/Lab4/diagrams/models.gif
[*] The more I think to it the less I see it feasible to have a unique equivalent viscosity. For instance, if the force is constant then each surface deformation is taking place following an exponential trend. Since the sum of two exponential functions does not give an exponential function, it seems you can't reduce the 2x{spring+dashpot} system to a {spring+dashpot} with equivalent properties.

I tried and it is not possible to obtain an equivalent viscosity in the case you are describing [1], except for a contact between two spheres that have the same material parameters (k,c). If the two spheres have the same parameters, the equivalent viscosity is cn = cn1cn2/(cn1+cn2) = c/2.
More genrally, this formulation is a good approximation if k2*delta_2 - k1*delta_1 negligible wrt the other force terms and delta_1 ~ delta_2

I understand the damping should not be zero if one particle of the two is purely elastic, however the actual formulation gives an arbitrary value (cn =c), which is higher than if we have contact between two particles with the same damping (cn = c/2). And this does not seems appropriate. So I am not sure which formulation is the most appropriate and if it is good to let it like that.

 [1] http://www.sccs.swarthmore.edu/users/06/dluong1/e12/Lab4/diagrams/models.gif model B with Mt=0.

I agree that the current state is not satisfying, yet I don't know how to improve it in a very general way.
If you find a good idea, then you can change this.

Hi Raphaёl,

feel free to make changes for the case, when one of cn==0, but
please, do not break the logic, when cn>0 for both elements.

Thanks,

Anton

I thought again at the problem, and I think we should stay consistent with the formulation cn = cn1cn2/(cn1 + cn2) all along. Because, if we make an exception for cn1 or cn2 = 0, the behavior will not be continue in a way. Think for example at keeping cn2 at a certain value and reducing slowly cn1 until 0, you have a discontinuity in the behavior at 0 which is not consistent. It is just a matter of being consistent, so, I think it is good to change and apply also the formulation cn = cn1cn2/(cn1 + cn2) at zero , and if nobody disagree I am going to do it.
Raphael

I just realized there is the same formulation for the stiffness (and the rolling resistance), should I change that also (same logic as above) ? (I am not sure for the rolling resistance as I am not using it. )
Something like :
if (cn1>0) or (cn2>0){
  phys->cn = cn1*cn2/(cn1+cn2 );
                 }
else {
               phys->cn = 0;
       }
and the same for cs, kn, and ks (and maybe mR) seems ok to you ?

Raphael

I hope it is ok.

Please, add a function, which returns the necessary value to
escape code duplication. Thanks.

Anton

Hi all,
I made the change for kn, ks, cn and cs and didn't change the rest. I tested and it works as expected. (nothing changed except the behavior at l1 or l2 = 0 (l being cn, cs, kn or ks) . I didn't touch mR because I don't precisely know what it is for. Feel free to do it if you think it should also be done.
I commited it, hope it will work fine.
Raphael

2013/12/4 Raphaël Maurin <email address hidden>:
> I made the change for kn, ks, cn and cs and didn't change the rest. I tested and it works as expected. (nothing changed except the behavior at l1 or l2 = 0 (l being cn, cs, kn or ks) .

We have this situation always during facet-sphere interaction [1].
Mass of one of element is 0, so the kn1 or kn2 will always be 0
(as well as ks, cn and cs).

So, we have to keep the previous behavior.

Anton

[1] https://github.com/yade/trunk/blob/master/pkg/dem/ViscoelasticPM.cpp#L37

One possibility is to keep Raphael's change for the case massMultiply=False.
Overall, I don't get the point of defining stiffness and other things for a facet if this facet will have null mass, implying stiffness=others=0.
I don't know the details of this code, but it sounds strange.