The graphs of position, velocity and acceleration are wrong.

1 view (last 30 days)
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

Accepted Answer

Mathieu NOE
Mathieu NOE on 7 Sep 2021
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

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!