Lagrangian for a double spring pendulum connected through a rigid bar

[Mishal0488]

Hi Guys

I am looking for some guidance with regards to a Lagrangian problem I am trying to solve.
Please refer to the attached documents.

Please neglect all the forcing functions for the time being. I am currently just trying to simulate the problem using initial conditions only

I have written Matlab code to simulate the problem, however my results are not as expected. the "theta2" variable always has a 0 solution, in some cases imaginary numbers form part of the solution. using initial conditions on the theta2 variable generally causes no solution as the ODE solvers on Matlab develop tolerance errors.

Note that the holonomic constraint whereby (b+r) is the subject of the formula is developed on Matlab using the symbolic package.

I believe the source of the issues is due to the holonomic constraint.
Can anyone provide me with some guidance on how to handle the system? and/or any mistakes which I have made. Could the system be represented in a more elegant manner?

Apologies for any funny terminology or poor representation of aspects, I am a structural engineer, Lagrangian mechanics isn't something that forms part of our undergraduate education hence this has been self taught.

function xp = Statespace3(t,x,r,L,Cs,Cp,Fv1,Fv2,Fh,w,w2,w3,k,m1,m2)

A = (x(5)+r);
three = A*cos(x(1))*cos(x(3));
four = A*sin(x(1))*sin(x(3));

R1 = L^2*sin(x(3))^2;
R2 = -A^2*cos(x(1))^2*sin(x(3))^2;
R3 = -A^2*cos(x(3))^2*sin(x(1))^2;
R4 = -2*A*L*cos(x(3))^2*sin(x(1));
R5 = 2*A^2*cos(x(1))*cos(x(3))*sin(x(1))*sin(x(3));
R6 = 2*A*L*cos(x(1))*cos(x(3))*sin(x(3));

one = sqrt(R1+R2+R3+R4+R5+R6);
B = L*sin(x(3))-one+three+four;

three_dot = x(6)*cos(x(1))*cos(x(3))-A*x(2)*sin(x(1))*cos(x(3))-A*x(4)*cos(x(1))*sin(x(3));
four_dot = x(6)*sin(x(1))*sin(x(3))+A*x(2)*cos(x(1))*sin(x(3))+A*x(4)*sin(x(1))*cos(x(3));

R1_dot = 2*x(4)*sin(x(3))*cos(x(3))*L^2;
R2_dot = -2*x(6)*A*cos(x(1))^2*sin(x(3))^2+2*x(2)*A^2*cos(x(1))*sin(x(1))*sin(x(3))^2-2*x(4)*A^2*cos(x(1))^2*sin(x(3))*cos(x(3));
R3_dot = -2*x(6)*A*cos(x(3))^2*sin(x(1))^2+2*A^2*x(4)*cos(x(3))*sin(x(3))*sin(x(1))^2-2*A^2*cos(x(3))^2*x(2)*sin(x(1))*cos(x(1));
R4_dot = -2*x(6)*cos(x(3))^2*sin(x(1))^2+2*x(4)*A*cos(x(3))*sin(x(3))*sin(x(1))^2-2*x(2)*A*cos(x(3))^2*sin(x(1))*cos(x(1));
R5_dot = 4*x(6)*A*cos(x(1))*cos(x(3))*sin(x(1))*sin(x(3))-2*A^2*x(2)*sin(x(1))^2*cos(x(3))*sin(x(3))-2*A^2*cos(x(1))*x(4)*sin(x(3))^2+sin(x(1))+2*A^2*cos(x(1))*x(2)*cos(x(1))^2*sin(x(3))+2*A^2*cos(x(1))*sin(x(1))*x(4)*cos(x(3))^2;
R6_dot = 2*x(6)*L*cos(x(1))*cos(x(3))*sin(x(3))-2*A*L*x(2)*sin(x(1))*cos(x(3))*sin(x(3))-2*A*L*x(4)*cos(x(1))*sin(x(3))^2+2*A*L*x(4)*cos(x(1))*cos(x(3))^2;

X_dot = (R1_dot+R2_dot+R3_dot+R4_dot+R5_dot+R6_dot);

b_dot = L*x(4)*cos(x(3))+three_dot+four_dot-X_dot/(2*(one));

%theta1 coordinate
F_t1 = A*Fh*sin(w*t)*cos(x(1));
theta1_dotdot = (F_t1-2*m1*x(2)*x(6)*A-9.81*m1*sin(x(1))*A-Cp*x(2))/(m1*A^2);

%theta2 coordinate
F_t2 = B*Fh*sin(w*t)*cos(x(3));
theta2_dotdot = (F_t2-2*m2*x(4)*b_dot*B-9.81*m2*sin(x(3))*B-Cp*x(4))/(m2*B^2);

%a coordinate
Fa = Fh*sin(w*t)*sin(x(1))+Fv1*sin(w2*t);
a_dotdot = (Fa+m1*x(2)^2*A-k*x(5)+9.81*m1*cos(x(1))-Cs*x(6))/m1;

%b coordinate
b_dotdot = (m2*x(4)^2*B+k*(B-r)*9.81*m2*cos(x(3)))/m2;
%state vector
xp = [x(2);
      theta1_dotdot;
      x(4);
      theta2_dotdot
      x(5)
      a_dotdot;];

 

[anuttarasammyak]

For the case of single spring pendulum
L=\frac{m}{2}(\dot{r}^2+r^2\dot{\theta}^2)-mgr sin\theta -\frac{k}{2}(r-r_0)^2
\frac{\partial L}{\partial q}-\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}
for q=r,$\theta$ are
mr\dot{\theta}^2 - mg \sin\theta -k (r-r_0) - m \ddot{r}
-m\dot{(r^2\dot{\theta})}- mgr cos\theta
where $r_0$ is natural length of spring.
By method of Lagrangean multiplier the equations of motions are
mr_1\dot{\theta_1}^2 - mg \sin\theta_1 -k (r_1-r_0) - m \ddot{r_1}+\lambda(t)\frac{\partial C}{\partial r_1}=0
mr_2\dot{\theta_2}^2 - mg \sin\theta_2 -k (r_2-r_0) - m \ddot{r_2}+\lambda(t)\frac{\partial C}{\partial r_2}=0
-m\dot{(r_1^2\dot{\theta_1})}- mgr_1 \cos\theta_1 +\lambda(t)\frac{\partial C}{\partial \theta_1}=0
-m\dot{(r_2^2\dot{\theta_2})}- mgr_2 \cos\theta_2 +\lambda(t)\frac{\partial C}{\partial \theta_2}=0
where
C=r_1^2+r_2^2-2r_1r_2cos(\theta_2-\theta_1)+2r_2L \sin\theta_2 -2r_1L\sin\theta_1
These four equations together with
C=0
which means the two spring pendulums are connected with a bar of length L.
should give us a solution.

I hope this would help you to check your way.

 

[Mishal0488]

Thank you.

With regards to the formulation of the Lagrangian we have done a similar thing with the difference being that my r1 and r2 are (a+r) and (b+r) however I think both a suitable.

I don't understand your holonomic constraint, can you please explain further?

Can the problem also be solved by simply making one of the variables of the holonomic constraint the subject of the formula (lets say r1)? We can then also take the time derivative of r1, thus using substitution we can reduce the number of DOF down to 3 instead of the current 4.

 

[anuttarasammyak]

Connection by bar is expressed as
(r_2 \cos\theta_2-r_1 \cos\theta_1)^2+(r_2 sin\theta_2+L-r_1 \sin\theta_1)^2=L^2
This is transformable to C=0 in post #2.

I have just got to make equations and not gone to how to solve them. Numerical solution might be necessary.

 

[anuttarasammyak]

Can the problem also be solved by simply making one of the variables of the holonomic constraint the subject of the formula (lets say r1)?

Another approach is to take x, y coordinates of the center of the bar and angle of bar wrt horizon as variables. We reduce parameters from four to three with the constraints built in.

 

[Mishal0488]

@anuttarasammyak, Thank you for the assistance. I managed to get the system to solve on Matlab, I have not had a chance to check the legitimacy of the results as yet.

The only difference between your suggestion and what I have is the manner in which the spring length is represented. Your method developed cleaner equations though due to there being less variables to worry about. Still can't figure out the issue with what I initially did...