弹簧质量阻尼系统数值算例—计算理论值

clc; clear; close all
%% 输入矩阵
M = 2.5e3*[1 0 0; 0 1 0; 0 0 1];
K = 1800000*[4 -1 0; -1 2 -1; 0 -1 3];
a = 0.001; b = 0.001; 
C = a*M + b*K;
%% 求解
[n, n] = size(M); %自由度数量；
[V, D] = eig(K, M);%广义特征值；
[ome ind] = sort(diag(D));%对特征值进行排序；
V = V(:, ind);%对振型进行排序；
xishu = diag(V'*M*V); %求解广义质量；
zhenxing = V*inv(sqrt(diag(xishu)))/(2*pi);%关于质量正交的振型；关于质量正交意思就是每个振型对应等效质量为1；
zizhenpinlv = diag(sqrt(zhenxing'*K*zhenxing)); %自振圆频率；
zhenxing1 = zhenxing./(zhenxing(1, :));
%% 振型就是变量zhenxing 一个频率对应一列；自振频率是zizhenpinlv
%% 阻尼比
dengxiaogangdu = diag(zhenxing'*K*zhenxing);
dengxiaozuni = diag(zhenxing'*C*zhenxing);
for i = 1:n
    zunibi(i) = dengxiaozuni(i)/(2*1*zizhenpinlv(i))*2*pi;% 每个振型对应的阻尼比，公式：kesia = C/（2Mw）
end
fprintf('自振频率为：%f\n', zizhenpinlv);
fprintf('阻尼比是: %f\n', zunibi*100);