求解由质量、弹簧、阻尼器和库仑摩擦组成的系统

问题描述

考虑以下系统:

图 1 - 质量、弹簧、阻尼器和库仑摩擦力(图片由 ，我有一个模糊的代码，我不知道如何完成:

from scipy.integrate import odeint将 numpy 导入为 np米 = 1.0k = 2.0c = 0.1亩 = 0.3穆克 = 0.2克 = 9.8vf = 0.01def eq(X, t, Xi):Ff = k * (Xi[0] - X[0]) + c * (Xi[1] - X[1]) # - m * dydt如果 np.abs(X[1]) 

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

如果您能帮助我了解对该系统进行数值求解的最规范方法是什么，我将不胜感激.提前致谢.

解决方案

另一种方法是使用 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((步骤，))dX[0] = v0ddX = np.zeros((步骤，))Ff = np.zeros((steps,))# FS = np.zeros((steps,))dt = t1/(步骤 - 1)对于 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])如果不是 (np.abs(dX[ii]) 

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

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

Consider the system below:

        Fig.1 - Mass, spring, damper and Coulomb frction (image courtesy of Wikimedia).

with a dynamic equation of:

where Ff is the Amontons-Columb friction defined as:

and consequently, the no-slip condition is defined as

Following this example, I have a vague code in mind which I don't know how to complete:

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)

where Xi0 is a step function. My main issue is that when I want to define Ff it depends on dydt which is to be defined later in that scope!

I would appreciate if you could help me know what is the most canonical way to numerically solve this system. Thanks in advance.

解决方案

another approach is just to use a for loop and calculate steps sequentially:

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

the result is shown as green in below plot:

I also changed the vf to 0.001. The results are different from odeint for some reason!

这篇关于求解由质量、弹簧、阻尼器和库仑摩擦组成的系统的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！