0

I'm trying to plot a diagram of Coulomb damping mass-spring model, see image bellow. mi*N is the damping friction force. enter image description here

Becuase the damping force switches direction depending on the direction of speed vector, we have 2 different equations of motion x(t), I named them x_a(t) and x_b(t). If you look closely on the diagram above, the use of x(t) depends on the period, for the first period which is from 0 to pi/omega_n x_a(t) is used, for the second period which is from pi/omega_n to 2*pi/omega_n we use x_b(t) and so on.

I was thinking to use a piece of code which would go something like

t = i * (pi/omega_n)
for i = 1, 3, 5, 7,...
x_a(t)
for i = 2, 4, 6, 8,...
x_b(t)

But I have no clue how to implement this. I managed to define the x_a(t) part of the code and plot it with the undampled model, see bellow.

import numpy as np
import matplotlib.pyplot as plt

#constants
k = 2    #(N/m), spring coef
m = 0.04  #(kg), mass
x0 = -0.1   #(m), preload
mi = 0.3    #(), dry dynamic friction coef. ABS-ABS
N = 0.3     #(N), normal contact force

f_tr = mi * N / k   #friction force/pring coef - equivalent distance
omega_0 = np.sqrt(k/m)

#time
t = np.linspace(0,5,100)

#undamped model
x_undamp = x0*np.cos(omega_0*t)
dx_undamp = -omega_0*x0*np.sin(omega_0*t)

#damped model
x_damp = (x0+f_tr)*np.cos(omega_0*t)-f_tr
dx_damp = -omega_0*(x0+f_tr)*np.sin(omega_0*t)

#time to x=0
t0 = np.arccos(f_tr/(x0+f_tr))/omega_0
print t0

#plotting
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.suptitle('position on left, velocity on right')
ax1.plot(t, x_undamp,'r', label='x_undamp')
ax1.plot(t, x_damp, 'b', label = 'x_damp')
ax2.plot(t, dx_undamp, 'r', label = 'dx_undamp')
ax2.plot(t, dx_damp, 'b', label = 'dx_damp')

#grids, titles, legends, axis labels
ax1.grid()
ax2.grid()
ax1.set_title('Position vs time')
ax2.set_title('Velocity vs time')
ax1.legend()
ax2.legend()
ax1.set_xlabel('t(s)')
ax1.set_ylabel('x(m)')
ax2.set_xlabel('t(s)')
ax2.set_ylabel('dx(m/s)')

plt.show()
user2882635
  • 133
  • 2
  • 19

1 Answers1

1

I'm gonna put plot a diagram of Coulomb damping mass-spring model solution here till someone comes with a better one:

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

m = 0.2
k = 2.0
c = 0.1
mus = 0.3
muk = 0.2
g = 9.8
vf = 0.01
v0 = 0.0
t1 = 10
sign = lambda x: np.tanh(100*x)

def Xi(t):
  if t < 1 :
    return 0
  else:
    return 1

vXi = np.vectorize(Xi)

def eq(X, t, Xi):
  Ff = k * (Xi(t) - X[0])
  
  if np.abs(X[1]) < vf and np.abs(Ff) < mus * m * g :
    Ff = k * (Xi(t) - X[0])
  else:
    Ff = sign(X[1]) * muk * m * g
  
  d2x = (k * (Xi(t) - X[0]) - Ff) / m
  return [X[1], d2x]

t = np.linspace(0, t1, 1000)
X0 = [v0, 0]
sol = odeint(func = eq, y0 = X0, t = t, args = (Xi, ), mxstep = 50000, atol = 1e-5)

plt.plot(t, sol[:, 0], 'r-', label = 'Output (x(t))')
plt.plot(t, vXi(t), 'b-', label = 'Input (xi(t))')
plt.ylabel('values')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()
Shahriar Hossain
  • 119
  • 2
  • 13
  • I suppose the code above is from the link (https://stackoverflow.com/questions/56754421/solving-a-system-of-mass-spring-damper-and-coulomb-friction) with viscous damping taken out. This is pretty new to me, so I'm approaching it step-by-step. I suppose Xi(t) is the driving force, constant in this case. Could you explain what vf and v0 are? – user2882635 Nov 09 '20 at 11:36