0
$\begingroup$

I am trying to solve a double spring-mass system that is positioned vertically (I currently have the pendulum at an angle of 0 in my code below). I derived the system of ODE's required, and I get this free body diagram and system of equations:

enter image description here enter image description here

The plot of the solution, given boundary conditions: r1_0=20, r2_0=40, theta_0=0:

enter image description here

The full Python code:

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import display, clear_output,HTML


#pendulum
b1=500
g=9.81

#mass 1
m1=5
k1=4
b2=0

#mass 2
m2=4
k12=10
b3=0

#Shift the mass positions
R1=0
R2=0

#Initial conditions
theta_0=0
omega_0=0

r1_0=20
v_r_1_0=0

r2_0=40
v_r_2_0=0


S_0=[r1_0,v_r_1_0,r2_0,v_r_2_0,theta_0,omega_0]

#DiffEq function
def dSdt(t,S):
    r1,r1dot,r2,r2dot,theta,thetadot=S
    return [r1dot,
            (-k1*(r1)-b2*r1dot-k12*((r1-r2))+m1*g*np.cos(theta))/m1,
            r2dot,
            (k12*((r1-r2))-b3*r2dot+m2*g*np.cos(theta))/m2,
            thetadot,
            (-thetadot*b1-m1*g*(r1)*np.sin(theta)-m2*g*(r2)*np.sin(theta))/(m1*((r1)**2)+m2*((r2)**2))
           ]
#Solve
t=np.linspace(0,60,400)
sol=sp.integrate.odeint(dSdt,y0=S_0,t=t,tfirst=True)
sol2=sp.integrate.solve_ivp(dSdt,t_span=(0,max(t)),y0=S_0,t_eval=t)

fig, axs = plt.subplots(2)
fig.suptitle('Mass Positions')
axs[0].set_ylim(-50,50)
axs[0].plot(t, sol2.y[0],label='r1 [m]')
axs[0].plot(t, sol2.y[2],label='r2 [m]')
axs[1].set_xlabel('Time [sec]')
axs[0].set_ylabel('Distance from pivot [m]')
axs[1].set_ylabel('Theta [deg]')
axs[0].invert_yaxis()
axs[0].legend()

axs[1].plot(t,sol2.y[4]*(180/np.pi),label='theta [rad]')

Again note that I currently have the pendulum at an angle of 0 (vertical).

How can I change the relative equilibrium positions of these masses so they aren't crossing each other? I need to include the weight in the equations because it will affect the motion while the pendulum is swinging.

$\endgroup$
2
  • $\begingroup$ Are you talking about the equilibrium positions, or the starting positions? It looks like your starting positions are far from equilibrium positions. Note that because of the nature of this problem, as long as the pendulum is moving the masses will also be moving along the pendulum: there is no equilibrium until the entire system is at rest. $\endgroup$
    – TimWescott
    Commented Jan 29 at 3:52
  • $\begingroup$ I meant the equilibrium positions. My answer below shows the solution that I came to. $\endgroup$ Commented Jan 29 at 15:43

3 Answers 3

2
$\begingroup$

At equilibrium $\ddot{r}_2=0,\ \dot{r}_2=0$. So, from the second equation, $$r_2-r_1 = \frac{m_2 g \cos\theta}{K_{12}}$$ For positive values of $m, g, \cos\theta, K_{12}$, $r_2 > r_1$. i.e. mass #2 is below mass #1 under equilibrium (i.e. static) conditions.

For linear equations, you can't prevent them crossing each other if the initial values of position and velocity are not favourable. Mathematically, "crossed" conditions satisfy the differential equations.

You can "relax the spring $K_{12}$ so "chances" of the masses crossing each other reduce.

$\endgroup$
9
  • $\begingroup$ Ok, thank you for the information. I initially tried to use a shift to lengthen the springs by an amount R1 and R2 by replacing r1 and r2 in the first two equations with r1-R1 and r2-R2, respectively. This resulted in the amplitude increasing instead of the equilibrium position. $\endgroup$ Commented Jan 28 at 13:40
  • $\begingroup$ That appears to be the way to directly influence the equilibrium positions. If I understood your comment correctly, $R_1$ and $R_2$ will then directly influence the equilibrium position (even if the spring stiffness is $\infty$). $\endgroup$
    – AJN
    Commented Jan 28 at 13:43
  • $\begingroup$ You can wait a day or two before accepting an answer. Maybe another answer which is more suitable might be posted. $\endgroup$
    – AJN
    Commented Jan 28 at 13:44
  • 2
    $\begingroup$ By setting theta to zero you have reduced the problem to the classic two spring mass in series system, aka a harmonic absorber amongst other things. I suspect your problem is that you should be using one measurement for the extension of the springs, that is R1=0 R2=0 is not helpful. Do everything in terms of one variable, x $\endgroup$ Commented Jan 28 at 21:36
  • 1
    $\begingroup$ @AJN these ODE's are nonlinear. Look at the Python code, the sine term is left in the coefficient. It will work to model large-angle oscillations. Also I did determine how to shift the equilibrium position, by simple substitution. $\endgroup$ Commented Jan 29 at 14:01
0
$\begingroup$

I played with it more today, and found that the equilibrium positions can be shifted by amounts $R_1$ and $R_2$, by using this equation:

$r_2−r_1=m_2*g*cos(θ)/K12$

Substituting the shift constants $R_1$ and $R_2$:

$r_2-R_2-(r_1-R_1)=m_2*g*cos(θ)/K12$

So a direct substitution into the force balance equations of the mass will shift the masses a specified amount away from the equilibrium position. The torque balance equation for the pendulum is left the same.

With $R_1=R_2=0$, with overdamping and zero initial velocity,the masses settle out at values as shown on the plot:

enter image description here

They can both be moved up 10 meters by using $R_1=R_2=-10$ enter image description here

The relative positions of the masses can also be altered by changing $R_1$ and $R_2$, such as in this plot: enter image description here

$\endgroup$
0
$\begingroup$

It doesn't look like you have accounted for all terms in the equations. And in systems such as these you also need to consider the unstretched spring positions.

Let's assume that the positions of the masses when the springs are unstretched are $l1$ and $l2$ from the pivot. Let's also assume that the position of the masses relative to $l1$ and $l2$ are $r1(t)$ and $r2(t)$ (with the downward direction considered as positive).

I'm going to use the Lagrangian approach and Mathematica for deriving the equations.

The positions of the two masses are

p1 = {(Subscript[l, 1] + Subscript[r, 1][t]) Sin[\[Theta][
  t]], -(Subscript[l, 1] + Subscript[r, 1][t]) Cos[\[Theta][t]]};
p2 = {(Subscript[l, 2] + Subscript[r, 2][t]) Sin[\[Theta][
  t]], -(Subscript[l, 2] + Subscript[r, 2][t]) Cos[\[Theta][t]]};

$$\left\{\left(l_1+r_1(t)\right) \sin (\theta (t)),\left(-l_1-r_1(t)\right) \cos (\theta (t))\right\}$$ $$\left\{\left(l_2+r_2(t)\right) \sin (\theta (t)),\left(-l_2-r_2(t)\right) \cos (\theta (t))\right\}$$

Their velocities are

v1 = D[p1, t];
v2 = D[p2, t];

$$\left\{\left(l_1+r_1(t)\right) \theta '(t) \cos (\theta (t))+r_1'(t) \sin (\theta (t)),\left(-l_1-r_1(t)\right) \theta '(t) (-\sin (\theta (t)))-r_1'(t) \cos (\theta (t))\right\}$$ $$\left\{\left(l_2+r_2(t)\right) \theta '(t) \cos (\theta (t))+r_2'(t) \sin (\theta (t)),\left(-l_2-r_2(t)\right) \theta '(t) (-\sin (\theta (t)))-r_2'(t) \cos (\theta (t))\right\}$$

The kinetic energy of the system is

T = 1/2 Subscript[m, 1]  v1 . v1 + 1/2 Subscript[m, 2]  v2 . v2 //Simplify

$$\frac{1}{2} \left(2 l_1 m_1 r_1(t) \theta '(t)^2+2 l_2 m_2 r_2(t) \theta '(t)^2+l_1^2 m_1 \theta '(t)^2+l_2^2 m_2 \theta '(t)^2+m_1 r_1(t){}^2 \theta '(t)^2+m_2 r_2(t){}^2 \theta '(t)^2+m_1 r_1'(t){}^2+m_2 r_2'(t){}^2\right)$$

The potential energy is

U = 1/2 Subscript[k, 1]  Subscript[r, 1][t]^2 +
  1/2 Subscript[k, 2] (Subscript[r, 2][t] - Subscript[r, 1][t])^2 + 
  Subscript[m, 1]  g  p1[[2]] + Subscript[m, 2]  g  p2[[2]]

$$g m_1 \left(-l_1-r_1(t)\right) \cos (\theta (t))+g m_2 \left(-l_2-r_2(t)\right) \cos (\theta (t))+\frac{1}{2} k_1 r_1(t){}^2+\frac{1}{2} k_2 \left(r_2(t)-r_1(t)\right){}^2$$

The Lagrangian is $L = T - U$.

From which we can assemble the equations of motion

eqns = Simplify@Thread[Table[D[D[L, q'[t]], t] - D[L, q[t]], 
  {q, {\[Theta], Subscript[r, 1], Subscript[r, 2]}}] == 
  {0, -Subscript[b, 1] Subscript[r, 1]'[t], -Subscript[b, 2]  Subscript[r, 2]'[t]}] 

$$l_1 m_1 \left(g \sin (\theta (t))+2 r_1(t) \theta ''(t)+2 r_1'(t) \theta '(t)\right)+l_2 m_2 \left(g \sin (\theta (t))+2 r_2(t) \theta ''(t)+2 r_2'(t) \theta '(t)\right)+g m_1 r_1(t) \sin (\theta (t))+g m_2 r_2(t) \sin (\theta (t))+l_1^2 m_1 \theta ''(t)+l_2^2 m_2 \theta ''(t)+m_1 r_1(t){}^2 \theta ''(t)+m_2 r_2(t){}^2 \theta ''(t)+2 m_1 r_1(t) r_1'(t) \theta '(t)+2 m_2 r_2(t) r_2'(t) \theta '(t)=0$$

$$b_1 r_1'(t)+k_1 r_1(t)+k_2 \left(r_1(t)-r_2(t)\right)=m_1 \left(g \cos (\theta (t))+l_1 \theta '(t)^2+r_1(t) \theta '(t)^2-r_1''(t)\right)$$

$$m_2 \left(g \cos (\theta (t))+l_2 \theta '(t)^2+r_2(t) \theta '(t)^2-r_2''(t)\right)+k_2 \left(r_1(t)-r_2(t)\right)=b_2 r_2'(t)$$

Solving for the equilibrium values we get two positions. The first is when the masses are directly below the pivot and the other is when they are directly above the pivot.

eql = Solve[
eqns /. Derivative[_][_][t] :> 0, {\[Theta][t], 
 Subscript[r, 1][t], Subscript[r, 2][t]}][[1 ;; 2]] /. C[1] -> 0

$$\left\{\left\{r_1(t)\to \frac{g \left(m_1+m_2\right)}{k_1},r_2(t)\to \frac{g \left(k_2 m_1+k_1 m_2+k_2 m_2\right)}{k_1 k_2},\theta (t)\to 0\right\},\left\{r_1(t)\to \frac{-g m_1-g m_2}{k_1},r_2(t)\to \frac{-g k_2 m_1-g k_1 m_2-g k_2 m_2}{k_1 k_2},\theta (t)\to \pi \right\}\right\}$$

The equilibrium positions using the parameter values you have and some values values for $l_1$ and $l_2$.

pars = {g -> 9.81, 
 Subscript[m, 1] -> 5, Subscript[k, 1] -> 4, Subscript[b, 1] -> 500, 
 Subscript[m, 2] -> 4, Subscript[k, 2] -> 10, Subscript[b, 2] -> 0, 
 Subscript[l, 1] -> 50, Subscript[l, 2] -> 80};

eql /. pars

$$\left\{\left\{r_1(t)\to 22.0725,r_2(t)\to 25.9965,\theta (t)\to 0\right\},\\ \left\{r_1(t)\to -22.0725,r_2(t)\to -25.9965,\theta (t)\to \pi \right\}\right\}$$

Simulating the system from some initial position of the masses below the pivot we see that $\theta$ is $0$. Because of the large stiffness $b_1$ there is hardly any oscillations in $m_1$, but it eventually reaches the equilibrium value of $50+22.07$. And $m_2$ is very oscillatory because the stiffness is $0$ and it reaches the equilibrium value of $80+25.99$.

sols = NDSolveValue[Join[eqns /. pars, {Subscript[r, 1][0] == 20, 
  Subscript[r, 1]'[0] == 0, Subscript[r, 2][0] == 20, 
  Subscript[r, 2]'[0] == 0, \[Theta][0] == 0, \[Theta]'[0] == 0}], 
  {\[Theta][t], Subscript[l, 1] + Subscript[r, 1][t], 
  Subscript[l, 2] + Subscript[r, 2][t]} /. pars, {t, 0, 500}];
Plot[#, {t, 0, 500}, PlotRange -> All] & /@ %

enter image description here

Zooming into the plots show the oscillations more clearly.

Plot[#, {t, 0, 10}, PlotRange -> All] & /@ sols

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.