MATLAB Answers

bode plot from mass, stiffness and damping matrix

37 views (last 30 days)
TG
TG on 24 Dec 2018
Commented: TG on 28 Dec 2018
Hi, I have calculated M0, K0 and D0 matrices which I want to plot in a bode diagram. However, I plot it in another way which can be seen below:
K0 = [6000000 600000 0 0
600000 300000 0 0
0 0 364560 0
0 0 0 261510];
M0 = [10.0000 0.5400 0.5000 0.5000
0.5400 1.3436 0.2600 0.2800
0.5000 0.2600 1.0000 0
0.5000 0.2800 0 1.0000];
D0=[-242.3846 -121.3745 -120.0000 -101.6000
-121.3745 -65.3490 -62.4000 -56.8960
-120.0000 -62.4000 -240.0000 0
-101.6000 -56.8960 0 -203.2000];
%% Calculating the eigenfrequencies of the system
[ur,omegaSqr]=eig(K0,M0);
omega=sqrt(diag(omegaSqr));
ur1norm=ur(:,1)/norm(ur(:,1),Inf);
ur2norm=ur(:,2)/norm(ur(:,2),Inf);
eigenf=omega/2/pi
%% plotting bode diagrams
%FREQUENCIES=linspace(0,200,20000);
FREQUENCIES=[0:0.5:eigenf(1)-0.5 eigenf(1)-0.5:0.001:eigenf(1)+0.5 ...
eigenf(1)+0.5:0.5:eigenf(2)-0.5 eigenf(2)-0.5:0.001:eigenf(2)+0.5 ...
eigenf(2)+0.5:0.5:eigenf(3)-0.5 eigenf(3)-0.5:0.001:eigenf(3)+0.5 ...
eigenf(3)+0.5:0.5:eigenf(4)-0.5 eigenf(4)-0.5:0.001:eigenf(4)+0.5 ...
eigenf(4)+0.5:0.5:200]; % more accuracy around eigenfrequency
OMEGA=FREQUENCIES*2*pi;
H11=zeros(length(OMEGA), 1);
H12=zeros(length(OMEGA), 1);
H21=zeros(length(OMEGA), 1);
H22=zeros(length(OMEGA), 1);
for oCount=1:length(OMEGA)
FRF=inv(-(OMEGA(oCount)^2)*M0+K0+OMEGA(oCount)*D0);
H11(oCount)=(FRF(1,1));
H12(oCount)=(FRF(1,2));
H21(oCount)=(FRF(2,1));
H22(oCount)=(FRF(2,2));
end
figure(1)
subplot(2,2,1)
plot(OMEGA./(2*pi),abs(H11))
grid on;
ylabel('$|H_{11}|$', 'Interpreter', 'latex', 'FontSize',14);
xlabel('$\Omega$ [Hz]','Interpreter', 'latex', 'FontSize', 14)
vline(eigenf);
subplot(2,2,2)
plot(OMEGA./(2*pi),abs(H12))
grid on;
ylabel('$|H_{12}|$', 'Interpreter', 'latex', 'FontSize',14);
xlabel('$\Omega$ [Hz]','Interpreter', 'latex', 'FontSize', 14)
vline(eigenf);
subplot(2,2,3)
plot(OMEGA./(2*pi),abs(H21))
grid on;
ylabel('$|H_{21}|$', 'Interpreter', 'latex', 'FontSize',14);
xlabel('$\Omega$ [Hz]','Interpreter', 'latex', 'FontSize', 14)
vline(eigenf);
subplot(2,2,4)
plot(OMEGA./(2*pi),abs(H22))
grid on;
ylabel('$|H_{22}|$', 'Interpreter', 'latex', 'FontSize',14);
xlabel('$\Omega$ [Hz]','Interpreter', 'latex', 'FontSize', 14)
vline(eigenf);
However, the tops of the graph are not on the places of the eigenfrequencies, only for the biggest peak. When I remove the damping response in the line
FRF=inv(-(OMEGA(oCount)^2)*M0+K0+OMEGA(oCount)*D0);
so replace D0 by a 0, the peaks are exactly on the places of the eigenfrequency, but no damping response is visible of course. I wonder whether there is another method to calculate this graph or plot a bode diagram when I have my M0, K0 and D0 matrices, or that I have made a problem somewhere. Thank you in advance!

Accepted Answer

Aquatris
Aquatris on 26 Dec 2018
First of all, I think your M K and D matrices are not for a spring mass damper system. They do not have any negative elements in them (if a spring is attached to both 1st and 2nd mass, individual motion of 1st mass and 2nd mass will have opposite effect on each other, hence sign difference). I recommend you check your problem formulation.
Having said that, the easiest way is to form a state space equation for the system. Below is the code for your system;
% x = [q1 q2 q3 q4 q1dot q2dot q3dot q4dot]'
% qi = position of ith mass
% qidot = velocity of ith mass
% xdot = Ax + Bu
% y = Cx
A = [zeros(4) eye(4);
-M0\K0 -M0\D0]; % system matrix
B = [zeros(4);% input matrix
M0\eye(4)]; % input 1 is force at 1st mass
% input 2 is force at 2nd mass etc
C = [eye(8)]; % output matrix
% output 1-4 is mass 1-4 position output
% output 5-8 is mass 5-8 velocity output
sys = ss(A,B,C,0);
bode(sys(3,4)) % bode diagram from input 4 to output 3

  1 Comment

TG
TG on 28 Dec 2018
Thank you, this is indeed what I wanted and I will check the matrices again!

Sign in to comment.

More Answers (0)

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by