how to solve this equation of motion?

4 views (last 30 days)
AJ
AJ on 21 Jan 2020
Commented: darova on 23 Jan 2020
Hi
I appreciate if anyone can help me to underestand if I can solve this EOM like a normal 2DOF system which I use ODE45 for it.
Here is the paper:
3.JPG
1.JPG
which:
2.JPG
function yp = func_forced_5DOF(t,y,M,ws,m1,m2,r,kx,ky,ksi,fx,fy,fsi,fd1,fd2,J,j1,j2,L1,L2)
a0 = (1 - ((ws^2)*Ts*Tr*Sig)^2) + ((ws^2)*((Tr + Ts)^2));
a1 = (2*(1 - ((ws^2)*Ts*Tr*Sig))*ws*Ts*Tr*Sig) - (2*(ws*((Tr + Ts)*Tr)));
a2 = Tr^2 + ((ws^2) * (Ts^2) * (Tr^2) *(Sig^2));
Te1 = (K*(ws-y(9))) / ((a2*(y(9)^2))+(a1*y(9))+a0);
Te2 = (K*(ws-y(10))) / ((a2*(y(10)^2))+(a1*y(10))+a0);
yp(1,:) = y(6,:);
yp(2,:) = y(7,:);
yp(3,:) = y(8,:);
yp(4,:) = y(9,:);
yp(5,:) = y(10,:);
yp(6,:) = ((m1*r*( ((y(9,:)^2) * cos(y(4,:))) + (yp(9,:)*sin(y(4,:))) )) - ( m2*r( ((y(10,:)^2) * cos(y(5,:))) + (yp(10,:)*sin(y(5,:))) )) - (fx*y(6,:)) - (kx*y(1,:))) / M;
yp(7,:) = ((m1*r*( ((y(9,:)^2) * sin(y(4,:))) - (yp(9,:)*cos(y(4,:))) )) + ( m2*r( ((y(10,:)^2) * sin(y(5,:))) - (yp(10,:)*cos(y(5,:))) )) - (fy*y(7,:)) - (ky*y(2,:))) / M;
yp(8,:) = ((-1*m1*r*L1*( ((y(9,:)^2) * sin(y(4,:))) - (yp(9,:)*cos(y(4,:))) )) + ( m2*r*L2*( ((y(10,:)^2) * sin(y(5,:))) - (yp(10,:)*cos(y(5,:))) )) - (fsi*y(8,:)) - (ksi*y(3,:))) / J;
yp(9,:) = (Te1 - (m1*r*( ( yp(7,:) * cos(y(4,:))) - (yp(6,:)*sin(y(4,:)) ))) - ( L1*( (yp(8,:) * cos(y(4,:))))) - (fd1*y(9,:)) ) / j1;
yp(10,:) = (Te2 - (m2*r*( ( yp(7,:) * cos(y(5,:))) - (yp(6,:)*sin(y(5,:)) ))) + ( L2*( (yp(8,:) * cos(y(5,:))))) - (fd2*y(10,:)) ) / j2;

Answers (1)

darova
darova on 21 Jan 2020
Here is an idea:
function main
% define constants
m1 = ...
m2 = ...
function dy = myode(t,u)
x = u(1);
dx = u(2);
%% ...
phi2 = u(9);
dphi2 = u(10);
% coefficient matrix
% d2x d2y d2psi d2phi1 d2phi2
A = [M 0 0 -m1*r*sin(phi1) m2*r*sin(phi2)
0 M 0 m1*r*cos(phi1) m2*r*cos(phi2)
0 0 J m1*r*L1*cos(phi1) ...
% ...
];
% constant matrix
B = [-fx*dx - kx*x + m1*r*dphi1^2*cos(phi1) - m2*r*dphi2^2*cos(phi2)
%..
];
dy(1:5,1) = u(2:2:10); % velocities
dy(6:10,1) = A\B; % accelerations
end
end
  6 Comments
darova
darova on 23 Jan 2020
Those values are something like M or kx. You can't guess them
Imagine you have something like that:
123.png
You want to see what will hapen if you loose green ball.
What will be the initial conditions (Velocity and angle)? Depends on what you want

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!