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

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

代码如下

%多自由度系统仿真

clc

clear

closeall

%%物理参数

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'), gridon

subplot(222)

plot(t, y1u2)

title('y1-u2'), gridon

subplot(223)

plot(t, y2u1)

title('y2-u1'), gridon

subplot(224)

plot(t, y2u2)

title('y2-u2'), gridon

%%解多自由度问题

[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

上面代码中的MDOFSolve函数如下

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

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

%求特征频率

symsw

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;

forr = 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)));

forr = 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