Fixed step size and end-of-step event

[Simone Hendriks]

Hi,

I have two questions about the integration of the EOM:

a) The system under study implements a force F=-A*sign(v), where A is a constant and v is the speed of the point where the force acts. The actual implementation smoothes the sign() function to avoid integration issues. The problem of the latter is that it introduces an additional parameter (a "smoothing stiffness") which is not always easy to set, and it leads anyway to small step sizes when the stiffness parameter is large. CPODES already performs better for this particular system, but my question is, would it make sense to use the unsmoothed sign() function and fix the integration step size? I am not particularly concerned about the sign() transition, but I need to avoid the integrator to reduce the step size. For the most part, this system has v around 0, which leads to too small step sizes for most of the integration. How would you approach this?

b) Unrelated to the previous question, is there a way to call a function after an integration step has completed? I am aware of the auto-update state variables, but for a particular application it would be easier to do it so.

Thank you in advance!

Best regards,
Simone

[Michael Sherman]

Hi, Simone.

Apologies for the delayed reply.

(b) There is currently no provision in Simbody to call a function at step completion. This would be a nice addition to the TimeStepper class in case anyone would like to submit a Simbody PR!

(a) This is a great question. It's not clear to me whether a fixed step would be stable -- worth a try. You can of course improve the integration speed by relaxing the drift parameter, at the expense of stiction that doesn't completely stick. That's the most common fix but might not work for you. One possibility (depending on the specifics of your application) would be to try Clay Anderson's recently added "ExponentialSpring" contact model. That uses extra state variables to deal with stiction in a non-stiff way that integrators like. It is limited in its applicability, though -- see the documentation in the header here. (You would need to build Simbody from source to get that.) I don't have any experience with that myself so don't know how effective that might be. Would be good to hear if you give it a try, though, and I'm sure Clay would be thrilled to hear from you.

Regards,
Sherm

[Simone Hendriks]

Hi Sherm,

Thank you very much for the reply, and apologies for my late reply as well.

(b) I'm not yet confident about extending Simbody, but please if anyone reading interested in this feature wants to submit a PR, it would certainly come in handy.

(a) Clay's approach is indeed very interesting! While the class is not directly applicable in my case (F=-A*sign(v) is actually a torque between two bodies, and not related to any contact station), the theory behind could apply. If I understood correctly, the frictional term consists of two parts: a linear spring-damper (the term that actually avoids drift) and a pure damper, both saturated to mu*N. The rest position of the spring is updated when full sliding takes place, which depends on a specified transition velocity. Aside from the fact that this formulation can avoid drift, I am still unsure about the efficiency improvement here, because there is still a transition velocity that for some applications needs to be very small, but I guess the addition of a fictive spring allows to consider greater transition speeds, thus making the system less stiff? I'd like to give it a try as time allows and post the results, it will also be interesting to compare it against a smoothed approach.

Thank you again for the valuable insights!

Best regards,
Simone

[Simone Hendriks]

Hi,

I just wanted to post a summary of the results/thoughts I got using different methods, in case it might help anyone. I tried the following to model a Coulomb-like friction force (with a sign function discontinuity and depending on a reaction force):

a) Using an auto-update state variable to calculate the reaction forces and calculating the friction force after each integration step with sign().
b) Using an auto-update state variable to calculate the reaction forces and calculating the friction force in calcForce() with a smoothed sign() function.
c) Sampling the reaction forces at constant intervals using a PeriodicEventReporter and calculating the friction force after each integration step with sign().
d) Sampling the reaction forces at constant intervals using a PeriodicEventReporter and calculating the friction force in calcForce() with a smoothed sign() function.
e) Using, I would say, a similar method to Clay's approach for the tangential forces, that is, setting up an additional linear spring-damper whose setpoint is updated upon full sliding (determined by a transition speed) using an auto-update state variable, combined with a).

When aiming at getting similar performance, all methods yield very close results, at least for my particular system, where the reaction forces vary quite smoothly. The problem with b), c), d) and e) is that additional parameters are required, which is not suitable for my particular problem because other people less experienced with multibody simulation must run the program with different configurations, meaning that these additional parameters cannot be fixed. I am still not sure which approach makes more sense from a theoretical point of view, but for example I would say that a) makes more sense than b), because in b) the reaction forces are assumed constant during an integration step, whereas the speed varies within a step computation.

If drift is a concern, I believe that using e) a set of parameters can be found that make the integration fast while avoiding drift, but again the (numerical) spring-damper parameters and the transition speed must be chosen carefully. It would certainly be interesting to find a way of computing these parameters depending on the system, but I was not capable of devising an strategy here.

Any additional thoughts would be very much appreciated!

Best regards,
Simone