有阻尼多自由度系统固有频率、阻尼的求取

对如下弹簧阻尼系统做的仿真，可以看到这是一个二阶系统，有阻尼存在，且阻尼是比例阻尼，可以用传统模态矩阵的方法求解，不过如果要让matlab画出多自由度系统的响应，需要使用状态方程辅助矩阵，具体详见代码。

http://s16/mw690/002XKryigy6Qj1Y84pN6f&690

    代码如下

% 多自由度系统仿真

clc

clear

close all

 

%% 物理参数

m = 200;

k = 980*1000;

c = 1.5*1000;

 

m1 = 1;

m2 = 1;

k1 = 1;

k2 = 1;

k3 = 0;

c1 = 1;

c2 = 1;

c3 = 0;

 

u1 = 1; % 激励

u2 = 1;

 

% 物理方程的矩阵

M = [m1 0;

    0 m2;]*m;

C = [c1+c2 -c2;

    -c2 c2;]*c;

K = [k1+k2 -k2;

    -k2 k2;]*k;

Z = zeros(size(M));

I = eye(size(M));

F = [u1; u2];

 

%% 状态方程辅助矩阵做仿真用

A = [C M;

    M Z];

B = [K Z;

    Z -M;];

G = [F;zeros(size(F))];

D = [Z I;

    -inv(M)*K -inv(M)*C;];

 

% 标准状态方程矩阵

Asys = D;

Bsys = inv(A);

Csys = [I Z];

Dsys = [Z Z];

sys = ss(Asys, Bsys, Csys, Dsys);

 

% 脉冲响应

ts = 0.001;

fs = 1/ts;

t = 0:ts:3;

y = impulse(sys, t);

 

y1u1 = y(:,1,1);

y1u2 = y(:,1,2);

y2u1 = y(:,2,1);

y2u2 = y(:,2,2);

figure

subplot(221)

plot(t, y1u1)

title('y1-u1'), grid on

subplot(222)

plot(t, y1u2)

title('y1-u2'), grid on

subplot(223)

plot(t, y2u1)

title('y2-u1'), grid on

subplot(224)

plot(t, y2u2)

title('y2-u2'), grid on

 

%% 解多自由度问题

[omega, phi, phin] = MDOFSolve(M, K);

Mr = simplify(phin'*M*phin);

Cr = simplify(phin'*C*phin);

Kr = simplify(phin'*K*phin);

 

fprintf('转化为正则空间的结果:\n');

fprintf('Mr=\n');

eval_r(Mr)

fprintf('Cr=\n');

eval_r(Cr)

fprintf('Kr=\n');

eval_r(Kr)

 

omegan = diag(Kr).^0.5;

zeta = diag(Cr)/2./omegan;

 

fprintf('固有频率(Hz)\n');

subs(omegan/2/pi)

fprintf('阻尼比\n');

subs(zeta)

 

    得到的结果如下

转化为正则空间的结果:

Mr=

ans =

                         1    -8.01360251500077e-017

    -8.01360251500077e-017                         1

Cr=

ans =

          2.86474508437579     1.65595689001804e-016

     1.65595689001804e-016          19.6352549156242

Kr=

ans =

          1871.63345512552     1.08189183481179e-013

     1.08189183481179e-013          12828.3665448745

固有频率(Hz)

ans =

          6.88542150158426

          18.0262675180169

阻尼比

ans =

        0.0331089636830301

        0.0866803922544587

http://s3/mw690/002XKryigy6Qj20qzHYe2&690

    上面代码中的MDOFSolve函数如下

% 多自由度模态矩阵的求取

function [omega, phi, phin] = MDOFSolve(M, K)

% 求特征频率

syms w

eq = det(K - w*M);

omega2 = solve(eq, w); % omega平方

[~,IX] = sort(subs(omega2));

omega2 = omega2(IX); % 排序

 

omega2 = real(eval_r(omega2));

 

% fprintf('特征频率的平方omega2:\n');

% pretty(omega2);

 

% 求模态矩阵

omega2u = unique(omega2);

[~,IX] = sort(subs(omega2u));

omega2u = omega2u(IX);

u = sym(zeros(size(M,1)));

ind = 1;

for r = 1:length(omega2u)

    ur = null(K-omega2u(r)*M); % 重特征值的情况，M为单位阵的话，直接满足要求；否则要另处理

    u(:, ind:ind+size(ur, 2)-1) = ur;

    

    ind = ind + size(ur, 2);

end

 

% fprintf('模态矩阵u:\n');

% pretty(u)

% fprintf('u''*M*u:\n');

% pretty(simplify(u'*M*u))

% fprintf('u''*K*u:\n');

% pretty(simplify(u'*K*u))

 

% 求正则模态矩阵

unorm = sym(zeros(size(M,1)));

for r = 1:size(u, 1)

    ur = u(:, r);

    Mr = ur'*M*ur;

    ur = ur/Mr^0.5;

    unorm(:, r) = ur;

end

 

% fprintf('正则模态矩阵unorm:\n');

% pretty(unorm)

% fprintf('unorm''*M*unorm:\n');

% pretty(simplify(unorm'*M*unorm))

% fprintf('unorm''*K*unorm:\n');

% pretty(simplify(unorm'*K*unorm))

 

% out

% omega = simplify(sqrt(omega2));

omega = sqrt(omega2);

phi = simplify(u);

phin = simplify(unorm);

end