Let’s assume you want to model a multibody system that includes friction. This seems like a quite reasonable idea to me and in this post, I’ll share what I’ve learned when doing exactly this in GNU Octave (or Matlab – doesn’t matter much since we’ll use quite basic functionality).

(small note: scroll down to the end of the blog post for a link to my simulation code)

## Friction in static and dynamic systems – a very, very short introduction

The coulomb equation provides a simplified formula for the friction forces. In the most simple variant, it is determined with a maximum transmittable torque or force at which sliding occurs.

- sticking force: F < Fmax
- sliding force: F = Fmax

The maximum Force Fmax can be calculated with the normal force F_N and the friction coefficient µ using the Formula Fmax = µ * F. So one can very easily calculate static systems (I consider it a standard topic in engineering mechanics, e.g. “calculate for a certain mechanism involving friction weather it is still in static equilibrium or it will start to move”). Now, let’s think about what happens in dynamics:

In analytical mechanics, the coulomb friction formula leads to piecewise linear equations (one equation for sliding and one for sticking condition – you just start with one of the equations and switch to the other one when the transition between sliding and sticking occurs; the state at the transition point is the new initial condition).

So things get a bit more complicated, but an otherwise linear system can still be solved analytically. However, I was always told that in numerics, friction leads to a lot of ugly problems. Most multibody software appears to use some kind of simplification – in many cases, sticking is approximated by some kind of damping behavior.

So what happens when you have a multibody system and try to implement friction? Let’s find out.

## Resources, initial considerations and our example system

I’ve spent some time trying different keyword combinations with Google and was not able to find a decent example on how to implement friction in a system of ODEs for numerical dynamics. The most helpful pages I found are these:

- Physical Modeling of Mechanical Friction in Simulink (Mathworks Newsletter, 05/2008)
- Clutch Modeling: Part 1 & Part 2 (blog “Will there be a whiteboard?”)

Both provide excellent overview on the problems we face and possible implementations of friction models for Simulink and Modelica. I however wanted to have a clutch model in Matlab / GNU Octave and was unable to find an example implementation. So I had to write my own code (what eventually led to this post).

As an example, we’ll consider a system of two masses which are connected with a clutch. Since the most obvious application example is a vehicle with a mechanical clutch, we’ll use rotational masses for our example. And for better comparisons and so on, we’ll more or less stick to the model structure from the “Clutch modeling: part 1 & 2”-example from above.

So this is the example structure we’ll try to model:

If we had a function to calculate the clutch torque (or force) as F_clutch ( angle_1, angle_2, vel_1, vel_2 ), the state function for our system would look like this:

F_clutch_ = F_clutch (phi_1, phi_2, vel_1, vel_2); F_mass1 = F_constant - F_clutch; F_mass2 = F_clutch - vel_2 * damping_parameter; xp = [ vel_1; vel_2; F_mass1 / inertia1 ; F_mass2 / inertia2 ];

Okay, this is not complete but you get the idea. The procedure is like this:

- Calculate forces.
- Output the change of the state vector for the solver.

It would be very, very simple. Now, of course we could implement the friction force as a function of the sliding velocity. This will not allow actual sticking but if we are okay with having a little slip in the sticking region, it should actually work quite well. To be honest, I never took the time to implement this – the main reason was that I wanted my model to perform reasonably well regarding calculation time – but both posts I’ve linked above show examples on how this can be done.

## Sticking and slipping in the simulation model

To get a system with sticking and slipping, we need to distinguish between these two operating modes. So our procedure now will look like this:

- Calulate all other forces.
- Check if the two masses have the same speed.

No? → Sliding occurs

Yes? Check if the torque between them is below the maximum transmittable torque to distinguish between slipping and sticking. - Depending on result from (2), either add the sliding forces to both masses if sliding occurs or treat them as one mass if sticking currently occurs.
- Output the change of the state vector for the solver.

To check the sliding/sticking torque in point 2, we can analytically solve the torque in the connection between the two masses. So the idea is like this: Assume a rigid connection and calculate the connection torque. If the calculated value is above the maximum clutch torque, calculate sliding instead.

My implemented code looks like this:

% calculating the friction torque (assuming a sticking state): T_friction = (T1*J2 - T2*J1) / (J1+J2); % Solving the friction condition between the two inertias. % A small tolerance will be added to cope with numerical inaccuracies in the variables. if or( T_friction>T_friction_max, phip1>phip2+1e-3 ) % current state: sliding into positive direction phipp1 = (T1 - T_friction_max) / J1; phipp2 = (T2 + T_friction_max) / J2; elseif or( T_friction<-T_friction_max, phip1<phip2-1e-3 ) % current state: sliding into negative direction phipp1 = (T1 + T_friction_max) / J1; phipp2 = (T2 - T_friction_max) / J2; else % current state: sticking phipp1 = (T1 + T2) / (J1 + J2); phipp2 = phipp1; % correction of speed differences due to the above error margin if 1 phip1=(phip1+phip2)/2; phip2=phip1; end end

Compared to the code example above (the forces-only variant), this is a lot messier. Note that I calculated the external Torques T1 and T2 earlier in my code before the part shown here. Two other things are also worth discussing.

First point: I’ve added a small error margin for the speed comparisons:

phip1 > phip2 + 1e-3 phip1 < phip2 - 1e-3

It appears that in numerics, exactness is unachievable anyway. If you remove the error tolerances, the solver will stop solving your equation (at least this is what happens to me). The folks in the Matlab newsletter mention the same thing: “If the slip becomes small enough (defining a locking velocity range near zero), the surfaces stick together.” Note that they define a “range near zero”, not exactly zero. So this appears to be one of the tradeoffs which result from the numerical implementation.

The second point is a direct result from this trandeoff: Since sticking occurs at speeds close to zero, there may still be a small speed difference between the two masses left – resulting in a drift during sticking. So I’ve added some code to compensate for this possible difference during sticking:

phip1=(phip1+phip2)/2; phip2=phip1;

Again, this is a tradeoff from numerical inaccuracies and we just have to live with it. At least with my implementation – maybe there are better options out there I haven’t found.

## Simulation results

The code shown above is implemented in the sample scenario described earlier. I’ve chosen pretty standard values for the inertias and torques on both masses:

% inertias (in kg*m^2) J1 = 1; J2 = 1; % external torques (in Nm): T1 = 20; T2 = -phip2*10; % maximum friction torque (in Nm): T_friction_max = 50*(t>1)*(t<2);

Basically, the model simulates 1 second with open clutch (no friction torque), then the clutch torque is closed and there should be first a synchronization phase and afterwards a sticking phase. At t = 2 s, the clutch is opened again.

This is exactly what happens with the results:

The code uses the ODE45-Solver (standard for Matlab, included in the odepkg-package from octave-forge for GNU Octave). In the lower part of the graph, I’ve plotted the step size from the solver. Note how the solver decreases the step size when a transition between sliding and sticking occurs. The good thing is that due to the error control mechanism, the solver automatically finds the transition point between sliding and sticking. On the other hand, this only works if you have set a pretty small error margin for your solver so if you try broader error settings, this may lead to aborted simulation or completely wrong results.

I’ll stop with this post here – it is long enough already and I believe I’ve discussed the most important points. To sum things up, here’s what I have learned:

- How to implement friction in a multibody ODE system. It’s a bit messy, but it works.
- Requirements for actual sticking/slipping transition in a model is that you can directly modify both speed and forces on the connected masses. So if you use a software that let’s you only modify torques or forces, you can not implement actual sticking.
- Another requirement is that the ODE solver matches to your model requirements. As an example, a solver with a constant stepsize does not really work with the code I’ve just presented (maybe there’s a better way?). With my implementation, the maximum solver error actually needs to be quite low to produce correct results.

## Further tasks and ideas

Here are some small tasks if you would like to do a bit of extra thinking on the topic:

- Implement a friction model without locking. Compare simulation time and result quality. What works better?
- How could a model be set up if several masses are connected with friction? Note that locking and sliding can happen within all connections of the model in various combinations.

## Model code

You can download my simulation model here. It runs in both Matlab and GNU Octave.

## 2 thoughts on “Modelling friction in 1d dynamics”