bode plot from mass, stiffness and damping matrix

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!

[madhan ravi]

bode() ?

[TG]

You can't use bode(K0,M0,D0)

Answer 1

[Aquatris]

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

[TG]

Thank you, this is indeed what I wanted and I will check the matrices again!