python - 求解质量、弹簧、阻尼器和库仑摩擦系统
问题描述
考虑以下系统:
图 1 - 质量、弹簧、阻尼器和库仑系数(图片由Wikimedia提供)。
有一个动态方程:
其中Ff
Amontons-Columb 摩擦定义为:
因此,no-slip
条件定义为
按照这个例子,我脑子里有一个模糊的代码,我不知道如何完成:
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
稍后在该范围内定义哪个!
如果您能帮助我了解数值求解该系统的最规范方法,我将不胜感激。提前致谢。
解决方案
另一种方法是使用 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
因某种原因而不同!
推荐阅读
- javascript - React JS Firebase:调用函数刷新数据时的奇怪行为
- amazon-qldb - 在 QLDB 中查找已删除数据的历史记录
- python - 如何在列操作中找到字符串模式?
- php - 未定义变量:提交表单时的 conn
- r - 在平衡数据框中创建一个变量,该变量是另一列的转置
- php - 修复 CSV 标头的数组键
- angular - Angular 8:PatchValue 不能与 ChangeDetector 和 UpdateValue 一起使用
- machine-learning - Rapid Miner 中的 k-means 质心图实际上是什么意思(用外行的话来说)?
- nsis - 模拟标准 Windows 弹出消息 哔声
- c++ - boost 信号和槽不在不同的线程中工作(使用 boost::asio::io_service)