模拟扭转弹簧

Question

可能比编程问题更多的物理问题... 基本上，我试图模拟膝盖扭转弹簧的机器人（请参阅图片），并测量将哪些力施加到膝关节上以及弹簧变形的多少。

机器人的图片

我编写的程序制作的图表具有正确的形状，但是对于某些人来说，弹簧角度似乎与我的时间段的大小相扩展的原因，这绝对是不正确的。

我已经在下面附上了我的代码，但这也是我认为我在做的事情的详细说明：

我认为除了引力外，只有一支力量施加到机器人上，并且在开始时使用该力持续时间，直截了当。从那力我计算关节的角速度。

之后，我进入一个循环，在每个时间步长我更新弹簧的力和角度。我计算了取决于角度的重力。我认为我的模型类似于倒置的摆，腿质量的术语可以忽略不计。然后，我计算弹簧施加的反向，并从两者之间的差异中获取总力。

我计算力量的速度变化，并更新速度。从速度开始，我也可以更改角度并更新角度。我重复该循环，直到达到时间限制为止。

其余代码仅用于绘制。

由于某种原因，角度的尺度随着我选择的时间段的大小而大大变化。但是，除此之外，这些图形看起来如我所期望的：短暂振荡然后以一个值安定。我看到的一个问题是，弹簧的力直接取决于时间步，而重力不取决于时间步。由于其他一切都取决于两种力之间的差异，因此这会改变量表。但是，我不知道如何解决此问题。

感谢您的帮助！

import math
from scipy import constants
import matplotlib.pyplot as plt


class Spring:
    def __init__(self, k, l=0.17, m=6.0):
        self.k = k
        self.l = l
        self.m = m
        self.alphaplt = []
        self.forceplt = []

    def impulse(self, force, duration, at_angle, stept):
        time = 0
        alpha = at_angle
        i = self.m * self.l ** 2
        vk = (force / self.m * duration) / self.l
        while time <= 100:
            tm = 0.5 * self.m * constants.g * self.l * math.sin(alpha)
            tf = self.k * vk * stept / self.l
            tg = tm - tf
            vk = vk + (tg * stept * self.l) / i
            alpha = alpha + vk * stept
            time += stept
            self.alphaplt.append(alpha)
            self.forceplt.append(tg)

        plt.plot(self.alphaplt)
        plt.show()
        plt.plot(self.forceplt)
        plt.show()


if __name__ == "__main__":
    s = Spring(0.75)
    s.impulse(90, 0.5, 0, 0.01)

Answer 1

如果您想查看此系统的动画模拟，请转到
在线vpython的网页，注册（免费），粘贴我在下面写的脚本并运行它。代码的一些解释以脚本中的注释包含在脚本中

Web VPython 3.2

# system of differential equations of the mathematical model of the physical system:
def f(state, param, F):
    m, L, k, theta_0, g = param
    a = state.x
    p = state.y
    s_a = sin(a)
    c_a = cos(a)
    da_dt = p / (s_a**2)
    dp_dt = c_a * p**2 / (s_a**3) - (k/(2*m*L**2)) * (2*a + theta_0 - pi) + ((g + F/m) / (2*L)) * s_a 
    return vector(da_dt, dp_dt, 0.) 

# one Runge-Kutta time-step propagation of the system's dynamics:
def propagate_RK4(f, state, param, F, dt):
    k1 = f(state, param, F)
    k2 = f(state + dt*k1/2., param, F)
    k3 = f(state + dt*k2/2., param, F)
    k4 = f(state + dt*k3, param, F)
    return state + dt*(k1 + 2*k2 + 2*k3 + k4) / 6.

# transforming the dynamical variable angle into the actual, physical configuration of the system 
def configure(struct, state):
    struct[0].pos = vector( L*sin(state.x), L*cos(state.x), 0. )
    struct[1].pos = vector( 0.,  2*L*cos(state.x), 0. )
    struct[2].axis = struct[0].pos 
    struct[3].pos = struct[0].pos
    struct[3].axis = struct[1].pos - struct[0].pos
    return None


# initialization and animation of the system's dynamics:

# constant parameters of the model: 
m=6.0
L=0.17
k=55.75
theta_0 = 2*pi/3
g = 9.8
F = 150.
param = (m, L, k, theta_0, g)
a_0 = pi/10

# fixed location of the lowest, ground level, joint: 
sphere(pos = vector(0., 0., 0.), radius = .005, color=color.red)
 
# initial positions of the moving middle joint, where the torsion spring is
joint_1 = sphere(pos = vector(L*sin(a_0), L*cos(a_0), 0.), radius = .01, color=color.red, make_trail=True)

# initial positions of the upper joint, where the main load (box) is positioned
joint_2 = sphere(pos = vector(0., 2*L*cos(a_0), 0. ), radius = .01, color=color.red, make_trail=True)

# initial configuration of the bar between ground joint and joint 1
bar_01 = cylinder(pos=vector(0,0,0), axis=vector(0, L,0), radius=0.003, color=color.orange)

# initial configuration of the bar between ground joint and joint 2
bar_12 = cylinder(pos=vector(0,L,0), axis=vector(0,L,0), radius=0.003, color=color.orange)

# initial state of the system, given by an initial angle a_0:
state = vector(a_0, 0., 0.)
t = 0.
dt = 0.01
force_duration = 0.7
frame_rate = 15
force_duration = force_duration / frame_rate

# drawing in yellow the locations of joint 1 and 2 at time t=0 
sphere(pos = vector(L*sin(a_0), L*cos(a_0), 0.), radius = .005, color=color.yellow)
sphere(pos = vector(0., 2*L*cos(a_0), 0. ), radius = .005, color=color.yellow)

# animation of the dynamics of the system while the constant external force is acting: 
while t <= force_duration:
    rate(frame_rate)
    state = propagate_RK4(f, state, param, F, dt)
    configure([joint_1, joint_2, bar_01, bar_12], state)
    t=t+dt

# drawing in yellow the locations of joint 1 and 2 at the time moment when the external force stops acting
sphere(pos = vector(L*sin(state.x), L*cos(state.x), 0.), radius = .005, color=color.yellow)
sphere(pos = vector(0., 2*L*cos(state.x), 0. ), radius = .005, color=color.yellow)

# animation of the dynamics of the system  after the constant external force has stopped acting: 
while force_duration < t <= 10:
    rate(frame_rate)
    state = propagate_RK4(f, state, param, 0, dt)
    configure([joint_1, joint_2, bar_01, bar_12], state)
    t=t+dt