Equation of motion mechanics

[a.mlw.walker]

Hi I'm using lagrange mechanics to try and come up with the equation of motion for the system in the image attached.

Two pistons attached to a fly wheel, and are set 90 degrees apart on the wheel. 1 piston is pushed with a force F(t).
The flywheel has inertia I.

After using the lagrange equation I have the three left hand sides of my equaitons of motion.

$$M_{1}\ddot{x_{1}}+C_{1}\dot{x}$$=

$$Ir^{2}\ddot{\Theta}+0r\dot{\Theta}$$=

$$M_{2}\ddot{x_{2}}+C_{2}\dot{x}$$=

But I am struggling to work out the forces for the right hand sides. The is the force pushing piston 1, the flywheel has a moment of inertia, (and a momentum?) and the second piston is well I'm not sure what force this one is equal to... please can someone give me some help here

 

[John_V22]

There's only one unknown in this problem, so you should only end up with one equation. Lagrange's equations uses energy balance of the kinetic and potential energy, with derivatives of the energy with respect to the fundamental variable. In this case, use theta as the unknown. Write x1 and x2 as functions of theta, x1 = L cos(theta), x2 = L cos(theta+90). Write the total kinetic energy (T) in terms of theta and the total potential energy (V=zero).

Then take derivatives of (T-V) with respect to theta dot, and the second derivative with respect to time. Also take the derivate with respect to theta. The derivative with respect to theta should be zero since there is no change in potential energy. This would produce the left hand side.

THe right hand side is the generalized force, F * dr/d theta = F * d x1 / d theta.

This should produce 1 equation, as a function of theta only.

 

[a.mlw.walker]

Thank you, I will try that

 

[a.mlw.walker]

I tried to solve it myself, but if x1 and x2 are the relative displacement of the two pistons respectively, how does Lcos(theta) represent that?

 

[John_V22]

Since derivatives with respect to time are needed, it may be a good idea to keep track of t

For M1
Let x1 be the distance from the circle center to M1. The triangle problem is side(L), side(R), angle (theta) type, and the unknown is the other side. Let gamma be the angle opposite of side L, and beta be the angle opposite of side R. This is solve by
1) use The Law of Sines first to calculate one of the other two angles;
sin(theta(t)) / L = sin(gamma) / R
gamma = asin(sin(theta(t)) / L * R)
2) then use the three angles add to 180° to find the other angle;
gamma + beta + theta(t) = 180
beta = 180 - theta(t) - gamma = 180 - theta(t) - asin(sin(theta(t)) / L * R)
3) finally use The Law of Sines again to find the unknown side, y1.
sin(beta) / y1 = sin(theta(t)) / L
x1 = sin(beta) * L / sin(theta) = sin(180 - theta - asin(sin(theta(t)) / L * R)) * L / sin(theta(t))

Here's a MATLAB script to determine the two lengths
% Equations of motion

% A
% B
% M1 O M2
%
% O, A, and B are points on a circular disk with moment of inertia I
% O is the center
% A is connectected to horizontal moving mass M1 with link of length L
% B is connectected to horizontal moving mass M2 with link of length L
% Angle M1 O A is theta(t)
% Angle AOB is a right triangle
% A horizontal force, F(t) is applied to M1

% What is the equation of motion

clc
close all;

syms t R L
theta = sym('theta(t)');
M1_O_A = theta;

% Find the horizontal distance from O to M1.

% Look at triangle M1_O_A
% Angle M1_O_A is theta(t)
% Length A M1 is L

% 1) use The Law of Sines first to calculate
% one of the other two angles;
% sin(M1_O_A) / L = sin(A_M1_O) / R

A_M1_O = asin(sin(M1_O_A) / L * R);

% 2) then use the three angles add to 180° to find the other angle;
% A_M1_O + O_A_M1 + M1_O_A = pi
O_A_M1 = pi - (A_M1_O + M1_O_A);

% 3) finally use The Law of Sines again to find the unknown side, x1.
% sin(O_A_M1) / x1 = sin(M1_O_A) / L
x1 =sin(O_A_M1) / (sin(M1_O_A) / L);

% Find the horizontal distance from O to M2.

% Look at triangle M2_B_O
% Angle B_O_M2
% B_O_M2 + M1_O_A + 90 = 180

B_O_M2 = pi - (M1_O_A + pi/2)
% Length B M2 is L

% 1) use The Law of Sines first to calculate
% one of the other two angles;
% sin(B_O_M2) / L = sin(O_M2_B) / R

O_M2_B = asin(sin(B_O_M2) / L * R);

% 2) then use the three angles add to 180° to find the other angle;
% B_O_M2 + O_M2_B + M2_B_O = pi
M2_B_O = pi - (B_O_M2 + O_M2_B);

% 3) finally use The Law of Sines again to find the unknown side, x2.
% sin(M2_B_O) / x2 = sin(B_O_M2) / L
x2 =sin(M2_B_O) / (sin(B_O_M2) / L);


x1 =
-(L*sin(- asin(1/L*R*sin(theta(t))) - theta(t)))/sin(theta(t))

x2 =
(L*sin(pi/2 - asin(1/L*R*sin(1/2*pi - theta(t))) + theta(t)))/sin(pi/2 - theta(t))


Use the expressions for theta, x1, and x2 to determine kinetic energy (T).
potential energy (V) is zero

Lagrangian is T - V

 

[a.mlw.walker]

Thanks John for the hasty reply, internet on the blink here making replying very difficult.

I have never used MATLAB to actually find an equation of motion before, but from reading your code, I'm not sure that's what you actually meant is it?

Firstly when I ran it, there was an error, I changed all the theta(t)'s in the x1 and x2 equations to theta's and it ran, giving me expressions for x1 and x2 in the command window.

Was that what was supposed to happen?

Secondly, I have only ever used the lagrangian by hand. When I calculated the EOM, I assumed that because L>>R, I could say that

x1 = rsin(theta)
x2 = rsin(theta + 90)

I agree it is an assumption, but much simpler to solve the lagrangian.

I have attached a pdf of my solving for the lagrangian EOM, i would be extremely grateful if you could check it to see whether I have done the correct thing.

I have rearranged the equation for theta_double_dot because this is how it is needed for ODE45 as I understand it.

By the way, the EOM i am using at the moment is this:

pdot(2) = (f+sin(p(1))*cos(p(1))*(m1*r^2-m2*r^2)*p(2)^2-C1*r^2*(cos(p(1)))^2-C2*r^2*(sin(p(1)))^2)/(m2*r^2*(sin(p(1)))^2+I*r^2+m1*r^2*(cos(p(1)))^2);

f is the force applied.
It is different from the one in the pdf I think, but it gives correct graphs except that the rotational velocity is 4e8, which is clearly wrong, I am expecting around 50 rad/s.

Thank you very much for your time

alex