MATLAB Answers

how to solve this equation of motion?

13 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;

  0 Comments

Sign in to comment.

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

Show 3 older comments
darova
darova on 22 Jan 2020
Your code looks ok. I tried this:
t = [0 0.5];
Z0 = [0 0 0 0 0 1 0 1 0 1]; % all angular velocities are ones
I also changed the sign before A(5,3) and coefficient A(5,5)
-m1*r*sin(u(7)) m1*r*cos(u(7)) -m1*r*L1*cos(u(7)) j1 0;...
-m2*r*sin(u(9)) m2*r*cos(u(9)) m2*r*L2*cos(u(9)) 0 j2];
And got this:
img2.png
AJ
AJ on 23 Jan 2020
Thanks a lot. I really appreciate.
Is there any way to guess that what range the initial conditions are?
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.

Sign in to answer this question.