bode plot from mass, stiffness and damping matrix

16 次查看（过去 30 天）
TG
TG 2018-12-24
评论： TG 2018-12-28
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!

请先登录，再回答此问题。

采纳的回答

Aquatris
Aquatris 2018-12-26
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Control System Toolbox 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by