MATLAB Answers

The graphs of position, velocity and acceleration are wrong.

조회 수: 6(최근 30일)
Ankit Joshi
Ankit Joshi 2021년 9월 3일
답변: Mathieu NOE 2021년 9월 7일
Hello Everyone. Following is the code I'm using to calculate the response of coupled MDOF system. But the graphs I am getting is undesirable as it should damp but instead it's overshooting. Please Help me. Thanks.
clear all; close all; clc
tspan=[0 2]; %time in seconds
ic=[0.1 0.2 0.1 -0.1 0 0 0 0];
[t,x]=ode45(@diff,tspan,ic);
xf1=x(:,1); xf2=x(:,2); xf3=x(:,3); xf4=x(:,4);
fx1=x(:,5); fx2=x(:,6); fx3=x(:,7); fx4=x(:,8);
m1=95; m2=945; m3=10; m4=80; %Mass Values in kgs
k1=730000; k2=50000; k3=1000; k4=10000; %Stiffness Values in kg/sec^2
c1=10; c2=8000; c3=1000; c4=750; %Damping Values in kg/m.sec^2
M=[m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4]; %Mass Matrix
K=[-(k1+k2) k2 0 0; k2 -(k2+k3) k3 0; 0 k3 -(k3+k4) k4; 0 0 k4 -k4]; %Stiffness Matrix
C=[-(c1+c2) c2 0 0; c2 -(c2+c3) c3 0; 0 c3 -(c3+c4) c4; 0 0 c4 -c4]; %Damping Matrix
A=[-(k1+k2)/m1 k2/m1 0 0; k2/m2 -(k2+k3)/m2 k3/m2 0; 0 k3/m3 -(k3+k4)/m3 k4/m3; 0 0 k4/m4 -k4/m4]; %Characteristic Matrix
[e,d]=eig(A); %Eigen Values and Vector
w=sqrt(abs(d)); %Natural Frequency Matrix
F=1320; %Input force in Newton
a1=(K(1,1)*x(:,1)+K(1,2)*x(:,2)+C(1,1)*x(:,5)+C(1,2)*x(:,6))/M(1,1); %Acc M1
a2=(K(2,1)*x(:,1)+K(2,2)*x(:,2)+K(2,3)*x(:,3)+C(2,1)*x(:,5)+C(2,2)*x(:,6)+C(2,3)*x(:,7))/M(2,2); %Acc M2
a3=(K(3,2)*x(:,2)+K(3,3)*x(:,3)+K(3,4)*x(:,4)+C(3,2)*x(:,6)+C(3,3)*x(:,7)+C(3,4)*x(:,8))/M(3,3); %Acc M3
a4=(K(4,3)*x(:,3)+K(4,4)*x(:,4)+C(4,3)*x(:,7)+C(4,4)*x(:,8))/M(4,4); %Acc M4
subplot(3,1,1)
hold on;
plot(t,xf1,'r')
plot(t,xf2,'g')
plot(t,xf3,'b')
plot(t,xf4,'black')
ylabel('Displacement (m)')
subplot(3,1,2)
hold on;
plot(t,fx1,'r')
plot(t,fx2,'g')
plot(t,fx3,'b')
plot(t,fx4,'black')
ylabel('Velocity (m/s)')
title('Disp, Vel & Acc wrt T')
subplot(3,1,3)
hold on;
plot(t,a1,'r')
plot(t,a2,'g')
plot(t,a3,'b')
plot(t,a4,'black')
legend('m1','m2','m3','m4')
ylabel('Acceleration (m/s^2)')
xlabel('Time (s)')
clc
disp('Natural frequencies (rad/sec)')
disp(w)
disp('Velocity (m/s)')
disp(max(x(:,5)))
disp(max(x(:,6)))
disp(max(x(:,7)))
disp(max(x(:,8)))
disp('Acceleration (m/s^2)')
disp(max(a1))
disp(max(a2))
disp(max(a3))
disp(max(a4))
function yt = diff(t,x)
m1=95; m2=945; m3=10; m4=80; %Mass Values in kgs
k1=730000; k2=50000; k3=1000; k4=10000; %Stiffness Values in kg/sec^2
c1=10; c2=8000; c3=1000; c4=750; %Damping Values in kg/m.sec^2
M=[m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4]; %Mass Matrix
K=[-(k1+k2) k2 0 0; k2 +(k2+k3) k3 0; 0 k3 -(k3+k4) k4; 0 0 k4 -k4]; %Stiffness Matrix
C=[-(c1+c2) c2 0 0; c2 -(c2+c3) c3 0; 0 c3 -(c3+c4) c4; 0 0 c4 -c4]; %Damping Matrix
A=[-(k1+k2)/m1 k2/m1 0 0; k2/m2 -(k2+k3)/m2 k3/m2 0; 0 k3/m3 -(k3+k4)/m3 k4/m3; 0 0 k4/m4 -k4/m4]; %Characteristic Matrix
[e,d]=eig(A) %Eigen Values and Vector
w=sqrt(abs(d)) %Natural Frequency Matrix
F=1320; %Input force in Newton
yt=zeros(8,1);
yt(1)=x(5); %Velocity M1
yt(2)=x(6); %Velocity M2
yt(3)=x(7); %Velocity M3
yt(4)=x(8); %Velocity M4
yt(5)=(K(1)*x(1)+K(2)*x(2)+C(1)*x(5)+C(2)*x(6))/M(1); %Acc M1
yt(6)=(K(5)*x(1)+K(6)*x(2)+K(7)*x(3)+C(5)*x(5)+C(6)*x(6)+C(7)*x(7))/M(6); %Acc M2
yt(7)=(K(10)*x(2)+K(11)*x(3)+K(12)*x(4)+C(10)*x(6)+C(11)*x(7)+C(12)*x(8))/M(11); %Acc M3
yt(8)=(K(15)*x(3)+K(16)*x(4)+C(15)*x(7)+C(16)*x(8))/M(16); %Acc M4
end
  댓글 수: 2
Ankit Joshi
Ankit Joshi 2021년 9월 6일
@darova I've attached the image with original question. Please find attached.

댓글을 달려면 로그인하십시오.

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 9월 7일
hello
I didn't find any mistake in your original equations , but I am lazy and I found your code quite long and complicated... why not make simply a state space model and use all available tools to run a time domain simulation (using either step , impulse or lsim ( lsim Simulate time response of dynamic systems to arbitrary inputs.))
this is my code that I double checked and seems to work as expected. You can easily add the eigenvalues computation of the A matrix if you want
also here I plotted only the acceleration responses of the 4 masses for each step input
good luck
clearvars; close all; clc
M1=95; M2=945; M3=10; M4=80; %Mass Values in kgs
K1=730000; K2=50000; K3=1000; K4=10000; %Stiffness Values in kg/sec^2
C1=10; C2=8000; C3=1000; C4=750; %Damping Values in kg/m.sec^2
[A,B,C,D] = model_4masses(M1,M2,M3,M4,C1,C2,C3,C4,K1,K2,K3,K4);
% 4 plots showing a step response on each mass (separetely)
% plotted data : acceleration
for ci = 1:4
figure(ci)
step(A,B,C,D,ci)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,B,C,D] = model_4masses(M1,M2,M3,M4,C1,C2,C3,C4,K1,K2,K3,K4)
% state space model of the stacked 4 masses - springs - dash pots
A = [0 1 0 0 0 0 0 0;
-(K1+K2)/M1 -(C1+C2)/M1 K2/M1 C2/M1 0 0 0 0 ;
0 0 0 1 0 0 0 0;
K2/M2 C2/M2 (-K3-K2)/M2 (-C3-C2)/M2 K3/M2 C3/M2 0 0 ;
0 0 0 0 0 1 0 0 ;
0 0 K3/M3 C3/M3 (-K4-K3)/M3 (-C4-C3)/M3 K4/M3 C4/M3 ;
0 0 0 0 0 0 0 1;
0 0 0 0 K4/M4 C4/M4 -K4/M4 -C4/M4];
B = [0 0 0 0 0 0 ;
1/M1 0 0 0 0 0 ;
0 0 0 0 0 0 ;
0 1/M2 0 0 0 0 ;
0 0 0 0 0 0 ;
0 0 1/M3 0 0 0 ;
0 0 0 0 0 0 ;
0 0 0 1/M4 0 0];
C = [-(K1+K2)/M1 -(C1+C2)/M1 K2/M1 C2/M1 0 0 0 0
K2/M2 C2/M2 (-K3-K2)/M2 (-C3-C2)/M2 K3/M2 C3/M2 0 0
0 0 K3/M3 C3/M3 (-K4-K3)/M3 (-C4-C3)/M3 K4/M3 C4/M3
0 0 0 0 K4/M4 C4/M4 -K4/M4 -C4/M4];
D = [1/M1 0 0 0 0 0 ;
0 1/M2 0 0 0 0 ;
0 0 1/M3 0 0 0 ;
0 0 0 1/M4 0 0];
end

추가 답변(0개)

Community Treasure Hunt

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

Start Hunting!

Translated by