求解质量、弹簧、阻尼器和库仑摩擦系统答案

【问题标题】：Solving a system of mass, spring, damper and Coulomb friction求解质量、弹簧、阻尼器和库仑摩擦系统
【发布时间】：2019-11-07 07:32:35
【问题描述】：

考虑以下系统：

图 1 - 质量、弹簧、阻尼器和库仑系数（图片由 Wikimedia 提供）。

有一个动态方程：

其中Ff 是阿蒙顿-哥伦布摩擦，定义为：

因此，no-slip 条件定义为

在this example之后，我脑子里有一个模糊的代码，我不知道如何完成：

from scipy.integrate import odeint
import numpy as np

m = 1.0
k = 2.0
c = 0.1
mus = 0.3
muk = 0.2
g = 9.8
vf = 0.01

def eq(X, t, Xi):
  Ff = k * (Xi[0] - X[0]) + c * (Xi[1] - X[1]) # - m * dydt

  if np.abs(X[1]) < vf and np.abs(Ff) < mus * m * g :
    Ff = k * (Xi[0] - X[0]) + c * (Xi[1] - X[1]) # - m * dydt
  else:
    Ff = -np.sign(X[1]) * muk * m * g
    pass

  dxdt = X[1]
  dydt = (k * (Xi[0] - X[0]) + c * (Xi[1] - X[1]) - Ff) / m
  return [dxdt, dydt]

t = np.linspace(0, 10, 1000)
Xi0 = np.piecewise(t, [t < 1, t >= 1], [0, 1])
X0 = [0, 0]
sol = odeint(eq, X0, t)

其中Xi0 是一个阶跃函数。我的主要问题是，当我想定义Ff 时，它取决于稍后将在该范围内定义的dydt！

如果您能帮助我了解数值求解该系统的最规范方法，我将不胜感激。提前致谢。

【问题讨论】：

我不确定F_f 的给定表达式是否正确。你有这个的来源吗？通常在静态系统中考虑固体摩擦。这实际上是一个关于物理的问题，而不是 python 代码......我会说只有在滑点的情况下才需要求解动态方程，否则方程是 dX/dt=(acceleration=0, velocity=0)...
我想我已经在下面的帖子中正确地实现了它。这篇文章中的代码有一些问题。
@Foad 如果以下答案中的实现优于问题中的实现，请编辑问题以包含它们（或它们的较短版本等）

标签： python numpy numerical-methods ode

【解决方案1】：

我会在这里放一个简化/临时的解决方案，直到有人提出更好的解决方案：

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()

结果是：

我使用this、this 和this 帖子来编写代码。我忽略了阻尼和惯性以简化问题。

【讨论】：

【解决方案2】：

另一种方法是使用for循环并按顺序计算步骤：

Y = np.piecewise(t, [t < t2, t >= t2], [0, 1])
dY = np.insert(np.diff(Y) / np.diff(t), 0 , v0, axis = 0)
X = np.zeros((steps,))
dX = np.zeros((steps,))
dX[0] = v0
ddX = np.zeros((steps,))
Ff = np.zeros((steps,))
# FS = np.zeros((steps,))
dt = t1 / (steps - 1)

for ii in range(1, steps):
  X[ii] = X[ii - 1] + dt * dX[ii - 1]
  dX[ii] = dX[ii - 1] + dt * ddX[ii - 1]
  Ff[ii] = k * (Y[ii] - X[ii]) #+ c * (dY[ii] - dX[ii])
  if not (np.abs(dX[ii]) < vf and np.abs(Ff[ii]) < mus * m * g) :
    Ff[ii] = np.sign(dX[ii]) * muk * m * g
  ddX[ii] = (k * (Y[ii] - X[ii]) - Ff[ii]) / m

结果在下图中显示为绿色：

我还将vf 更改为0.001。由于某种原因，结果与odeint 不同！

【讨论】：

【解决方案3】：

写出这样一个系统的方程并不明显。而且解决起来也不容易。

如果可以解除Python约束，我建议使用OpenModelica来解决这个问题。在modelica组件库中，你有元素

.Modelica.Mechanics.Translational.Components.MassWithStopAndFriction

可用于模拟干摩擦。

【讨论】：