Simulating Torsion Springs

Question

First of all, I'm just a measly IT student and have very limited physics knowledge, so sorry, if my mistake is super obvious.

So basically I am trying to simulate a robot with a torsion spring in its knee (see picture) and measure what forces are applied to the knee joint and how much the spring deforms.

The program I've written produces graphs that have the correct shape, but for some reason the angle of the spring seems to be scaling with the size of my timesteps, which is definitely not correct.

I have attached my code below, but here is also a detailed explanation of what I think I am doing:

I assume that besides gravitation, there is only one force ever applied to the robot, and that force is applied at the start for a certain duration and is aimed straight down. From that force I calculate the angular velocity of my joint.

After that, I enter a loop where in every timestep I update the forces and the angle of the spring. I calculate the gravitational force dependent on the angle. I assume that my model resembles an inverted pendulum, and that the term for the mass of the leg is negligible. Then I calculate the counterforce exerted by the spring, and get the total force from the difference between the two.

I calculate the velocity change from the force, and update the velocity. From the velocity I get the change in angle and update the angle as well. I repeat that loop until a time limit is reached.

The rest of the code is just for plotting.

For some reason, the scale of the angle changes drastically with the size of the timesteps I choose. Other than that, however, the graphs look as I would expect: Oscillating briefly and then settling on one value. One problem I see is that the force of the spring directly depends on the timestep, while the gravitational force does not. Since everything else depends on the difference between the two forces, this changes the scale. However, I do not know how to fix this.

Thanks for any help!

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

The differential equation for this model should be something like

\begin{align} \frac{d}{dt} \left(\sin^2(\alpha) \frac{d\alpha}{dt}\right) \,=\, \sin(\alpha)\cos(\alpha)\left(\frac{d\alpha}{dt}\right)^2 -\, \frac{k}{2ml^2}\big(\,2\,\alpha + \theta_0 - \pi\,\big) \,+\, \left(\frac{g}{2l} + \frac{F}{2ml}\right)\sin(\alpha) \end{align}

which can be written as a first order system of differential equations

\begin{align} &\frac{d\alpha}{dt} \,=\, \frac{1}{\sin^2(\alpha)}\, p\\ &\frac{dp}{dt} \,=\, \frac{\cos(\alpha)}{\sin^3(\alpha)}\,p^2 \,-\, \frac{k}{2ml^2}\big(\,2\,\alpha + \theta_0 - \pi\,\big) \,+\, \left(\frac{g}{2l} + \frac{F}{2ml}\right)\sin(\alpha) \end{align}

The location of the joint (i.e. the rightmost point) has current coordinates \begin{align} &x\,=\, l\,\sin(\alpha)\\ &y\,=\, l\,\cos(\alpha) \end{align} The location of the bottom point on the box has current coordinates \begin{align} &x\,=\, 0\\ &y\,=\, 2\,l\,\cos(\alpha) \end{align}

The derivation of the equations of motion can be carried out by writing the Lagrangian of the system, then switching to Hamiltonian formulation.

Derivation. We start with writing the Lagrangian of the system $O \, \vec{i} \, \vec{j}$ in an inertial coordinate system attached with the ground. The origin of said system is at the ground-level joint, the x-axis is the horizontal axis aligned with the ground and the y-axis perpendicular to the ground. The Lagrangian of the system in an inertial coordinate system is the difference between the kinetic energies of all mass-points minus all potentials of the mass-points of the system. Here, only the top joint has mass assigned to it, so there is only one mass-point of mass $m$. It is not important to the dynamics of the model, but the location of the torsion spring joint is $$x\, \vec{i} \,+\, y\,\vec{j} \,=\, L\,\sin(\alpha) \, \vec{i} \,+\, L\,\cos(\alpha)\,\vec{j}$$ and the location of the top joint, the one with mass, is $$y\,\vec{j} \,=\,2\,l\,\cos(\alpha)\,\vec{j}$$ Velocity of the top joint: $$\frac{dy}{dt}\,\vec{j} \,=\,-\,2\,l\,\sin(\alpha)\,\frac{d\alpha}{dt} \,\,\vec{j}$$ In the inertial coordinate system, the kinetic energy and the potentials of the various forces are

kinetic energy: $$T \,=\,\frac{m}{2}\,\left(\frac{dy}{dt}\right)^2$$ torsion spring potential: $$U_{\text{spring}} \, =\, \frac{k}{2} \, \big(\,2\,\alpha + \theta_0 - \pi\,\big)^2$$ gravitational potential: $$U_{\text{gravity}} \, =\, mg\,y$$ finally the constant force applied vertically down to top joint is $$\vec{F}_{\text{external}}\,=\,F\,\vec{j}$$ and since this force is constant, it has a potential $$U_{\text{external}} \,=\, F\,y$$
Consequently, the Lagrangian in $\, O \, \vec{i} \, \vec{j} \,$ coordinates is simply $$\mathcal{L} \, =\, T \,-\, U_{\text{spring}} \, -\, U_{\text{gravity}} \,-\, U_{\text{external}}$$ $$\mathcal{L} \,=\, \frac{m}{2}\,\left(\frac{dy}{dt}\right)^2 -\, \frac{k}{2} \, \big(\,2\,\alpha + \theta_0 - \pi\,\big)^2 \, -\, \big(mg + F\big)\,y$$ Switching to (non-inertial, curvilinear) holonomic coordinate $y \,=\, 2 l\,\cos(\alpha)$ The Lagrangian becomes $$\mathcal{L} \,=\, \frac{m}{2}\,\left(2l\,\sin(\alpha)\frac{d\alpha}{dt}\right)^2 -\, \frac{k}{2} \, \big(\,2\,\alpha + \theta_0 - \pi\,\big)^2 \, -\, 2l\,\big(mg + F\big)\,\cos(\alpha)$$ $$\mathcal{L} \,=\, \frac{4ml^2}{2}\,\sin^2(\alpha)\,\left(\frac{d\alpha}{dt}\right)^2 -\, \frac{k}{2} \, \big(\,2\,\alpha + \theta_0 - \pi\,\big)^2 \, -\, 2ml\,\Big(\,g + \frac{F}{m}\,\Big)\,\cos(\alpha)$$ Consequently, the equation of motion of the system is the Euler-Lagrange equation derived from the Lagrangian: $$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\alpha}}\right) \,=\, \frac{\partial \mathcal{L}}{\partial {\alpha}}$$ which, written explicitly after dividing by some of the constants, is \begin{align} \frac{d}{dt} \left(\sin^2(\alpha) \frac{d\alpha}{dt}\right) \,=\, \sin(\alpha)\cos(\alpha)\left(\frac{d\alpha}{dt}\right)^2 -\, \frac{k}{2ml^2}\big(\,2\,\alpha + \theta_0 - \pi\,\big) \,+\, \left(\frac{g}{2l} + \frac{F}{2ml}\right)\sin(\alpha) \end{align} From here, set $$p \,=\, \sin^2(\alpha) \frac{d\alpha}{dt}$$ solve for $p$ and you get the Hamiltonian version of the dynamics

\begin{align} &\frac{d\alpha}{dt} \,=\, \frac{1}{\sin^2(\alpha)}\, p\\ &\frac{dp}{dt} \,=\, \frac{\cos(\alpha)}{\sin^3(\alpha)}\,p^2 \,-\, \frac{k}{2ml^2}\big(\,2\,\alpha + \theta_0 - \pi\,\big) \,+\, \left(\frac{g}{2l} + \frac{F}{2ml}\right)\sin(\alpha) \end{align}

This kind of system is easy to integrate numerically using some explicit numerical integrating method like Runge Kutta.