bode plot from mass, stiffness and damping matrix

34 visualizzazioni (ultimi 30 giorni)
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!

Risposta accettata

Aquatris
Aquatris il 26 Dic 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 Commento
TG
TG il 28 Dic 2018
Thank you, this is indeed what I wanted and I will check the matrices again!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Get Started with Control System Toolbox in Help Center e File Exchange

Prodotti


Release

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by