1
$\begingroup$

I have the equation of motion:

$a = \frac{1}{m} \left( F_1 - F_2 -F_{fric} \right)$

where a is the acceleration and m is the mass of an object. $F_{fric}$ is the friction calculated as: $F_{fric} = B \cdot v + F_{coulomb} \cdot sign(v)$

Where v is velocity and B is viscous friction coefficient.

I have a problem when implementing this in Simulink due to the coulomb friction as this increases the simulation time significantly and causes oscillations.

If I don't include the Coulomb term in the dynamic the position goes straight to the upper limit.

Can someone help me come across this problem?

To illustrate a picture of the Simulink structure and the scope of the forces:

enter image description here

enter image description here

Zoomed

enter image description here

Improve this question
$\endgroup$
4
  • $\begingroup$ Where did you get the formula $F_{coulomb} \cdot sign(v)$. This Wikipedia page doesn't seem to have this formula. As hinted by the existing answer, you need a better friction model. Wikipedia. $\endgroup$
    – AJN
    Commented Oct 27, 2022 at 17:13
  • $\begingroup$ This is the stribeck friktion model without the stribeck effect hindawi.com/journals/mpe/2012/432129/fig4 $\endgroup$
    – Nil
    Commented Oct 27, 2022 at 19:41
  • $\begingroup$ The attached paper clearly mentions below figure number 4, that, the signum function is used only for non zero values of velocity. They give a separate method to calculate friction force for the case where velocity is zero. Your implementation doesn't use this additional formula. $\endgroup$
    – AJN
    Commented Oct 28, 2022 at 2:10
  • $\begingroup$ Well, I only have the values for the dry friction and the viscous friction $\endgroup$
    – Nil
    Commented Oct 28, 2022 at 7:16

2 Answers 2

Reset to default
2
$\begingroup$

The sign function gives zero friction force for zero velocity, which is generally not correct (and leads to the artificial oscillations, as the integrator is made for smooth systems and cannot handle discontinuities), since for sticking the friction force is a constraint (absolute value between zero and maximum given by static friction), unless it exceeds the maximal value (static friction). Either you implement a boolean logic checking for transitions from sticking to slipping and vice versa (possibly available in SimMechanics?), or you use a continuous (and thus less accurate) regularization, e.g. an arctan- instead of a sign-function. Using MatLab instead of SimuLink you may apply event-detection, which is explained by the bouncing ball in the documentation (instead of collisions you would then have to check for stick-slip transitions).

Improve this answer
$\endgroup$
1
  • $\begingroup$ I am not sure I follow here ? how would you make this ? $\endgroup$
    – Nil
    Commented Oct 27, 2022 at 16:31
1
$\begingroup$

This is how it (time stepping, an alternative to event-detection algorithms) can be done in Python (I am afraid having neither SimuLink nor MatLab at hand, but the algorithm should become clear)

import numpy as np
import matplotlib.pyplot as plt

q0 = 0   # initial position
p0 = 0  # initial momentum
t_start = 0   # initial time
t_end = 10   # end time
N = 500 # time points
m = 1   # mass
k = 1   # spring stiffness

muN = 0.5   # friction force (slip and maximal stick)
omega = 1.5   # forcing radian frequency [RAD]
Fstat = 0.1   # static component of external force
Fdyn = 0.6   # amplitude of harmonic external force
F = lambda tt,qq,pp: Fstat + Fdyn*np.sin(omega*tt) - k*qq - muN*np.sign(pp)  # total force, note sign(0)=0 used to disable friction
zero_to_disable_friction = 0

omega0 = np.sqrt(k/m)
print("eigenfrequency   f = {} Hz;   eigen period   T = {} s".format(omega0/(2*np.pi), 2*np.pi/omega0))
print("forcing frequency   f = {} Hz;   forcing period   T = {} s".format(omega/(2*np.pi), 2*np.pi/omega))

time = np.linspace(t_start, t_end, N)   # time grid
h = time[1] - time[0]   # time step
q = np.zeros(N+1)   # position
p = np.zeros(N+1)   # momentum
absFfriction = np.zeros(N+1)

q[0] = q0   
p[0] = p0
for n, tn in enumerate(time):
   
    p1slide = p[n] + h*F(tn, q[n], p[n])   # end-time momentum, assuming sliding
    q1slide = q[n] + h*p1slide/m   # end-time position, assuming sliding
    
    if p[n]*p1slide > 0:   # sliding goes on
        q[n+1] = q1slide
        p[n+1] = p1slide
        absFfriction[n] = muN
        
    else:
        q1stick = q[n]   # assume p1 = 0 at t=tn+h
        Fstick = -p[n]/h - F(tn, q1stick, zero_to_disable_friction)    # friction force needed to stop at t=tn+h
        if np.abs(Fstick) <= muN:
            p[n+1] = 0   # sticking
            q[n+1] = q1stick
            absFfriction[n] = np.abs(Fstick)
        else:  # sliding starts or passes zero crossing of velocity
            q[n+1] = q1slide   # possible refinements (adapt to slip-start or zero crossing)
            p[n+1] = p1slide
            absFfriction[n] = muN
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.