新浪博客

[转载]Matlab 状态空间法求解阻尼系统固有频率

2015-01-15 13:13阅读：

http://blog.sina.cn/dpool/blog/u/2136851840

原文地址：Matlab 状态空间法求解阻尼系统固有频率
原文作者：长笛人倚楼Gloria

一、基本原理 1. 特征值
对于一个n自由度的阻尼系统，数学上可以以建立n个2阶常微分方程
MX'' + CX' + KX = F (1)
这时可以把位移和速度都看做独立的变量，将原来n个二阶常微分方程作为2n个一阶常微分方程处理。对于特征值问题，应有
(x = phi e^{lambda t}) (2)
将(2)式带入(1)式可以得到
((lambda^2M + lambda C + K)phi = 0) (3)
求解出来(lambda)就是其特征值。
2. 状态空间转换
现对式(1)添加如下方程
Mx' - Mx' = 0 (4)
(1)与(4)合写为
Ay' + By = G (5)
设
Z = zeros(size(M));
I = eye(size(M));
则有
A = [C M;
M Z];
B = [K Z;
Z -M;];
G = [F;zeros(size(F))];
y = [x'; x];


对于状态空间而言，应有X' = Asys X + Bsys U(t)(6)
将(6)式与(5)式比较可知，Asys = -inv(A)*B
Asys =
[Z I;
-inv(M)*K -inv(M)*C;];
我们知道，对于状态空间而言，Asys矩阵代表着系统的性质，因此，对Asys矩阵求特征值，就可以得到阻尼系统的特征值，进而得到固有频率和模态。
3. Asys特征与系统特征的关系
（1）A的特征向量与系统模态
设Asys的特征向量为psi，系统模态为phi
由于y = psi * e^(lambda*t);
x = phi * e^(lambda * t);
比较可知
psi = [lambda*phi; phi];
（2）A的特征值与系统固有频率
lambda的一般值是一对共轭复数
(lambda = -n pm iomega_d)
其中，n及(omega_d)即为系统作自由衰减震动的衰减系数与振动频率，若设
(omega_d = omega_i sqrt{1 - zeta^2}, n =zetaomega_i)
可知
(omega_i = sqrt{omega_d^2 + n^2}, zeta = frac{1}{sqrt{1 + (frac{omega_d}{n})^2}})
二、实例
例1：不考虑阻尼，进行计算
[转载]Matlab <wbr>状态空间法求解阻尼系统固有频率
利用振动力学方法人工计算结果为
[转载]Matlab <wbr>状态空间法求解阻尼系统固有频率
[转载]Matlab <wbr>状态空间法求解阻尼系统固有频率
下面以状态空间法，通过matlab求解，为简便，设k = 1， m = 1，代码如下：
M = [1, 0; 0, 2];
C = zeros(size(M));
K = [2, -1; -1, 3];

A = [zeros(size(M)), eye(size(M, 1)); -inv(M)*K, -inv(M)*C];
[V, D] = eig(A)
absD = abs(D)

结果如下：
V =
Columns 1 through 2
-0.0000 - 0.4781i -0.0000 + 0.4781i
0.0000 + 0.2390i 0.0000 - 0.2390i
0.7559 + 0.0000i 0.7559 + 0.0000i
-0.3780 + 0.0000i -0.3780 - 0.0000i
Columns 3 through 4
0.0000 + 0.5000i 0.0000 - 0.5000i
0.0000 + 0.5000i 0.0000 - 0.5000i
-0.5000 + 0.0000i -0.5000 + 0.0000i
-0.5000 + 0.0000i -0.5000 - 0.0000i
D =
Columns 1 through 2
-0.0000 + 1.5811i 0.0000 + 0.0000i
0.0000 + 0.0000i -0.0000 - 1.5811i
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i
Columns 3 through 4
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 1.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 - 1.0000i
absD =
1.5811 0 0 0
0 1.5811 0 0
0 0 1.0000 0
0 0 0 1.0000
可以看出，与上面手工计算的结果是一致的（这里的特征向量标准化过）。
例2：考虑阻尼，进行计算。
这个例子前面的博文有提到过http://blog.sina.cn/dpool/blog/s/blog_84024a4a0101e9lh.html?vt=4，计算结果为
固有频率(Hz)
ans =
6.88542150158426
18.0262675180169
[转载]Matlab <wbr>状态空间法求解阻尼系统固有频率
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;

% 物理方程的矩阵
M = [m1 0;
0 m2;]*m;
C = [c1+c2 -c2;
-c2 c2;]*c;
K = [k1+k2 -k2;
-k2 k2;]*k;

A = [zeros(size(M)), eye(size(M, 1)); -inv(M)*K, -inv(M)*C];
[V, D] = eig(A)
abs(D / 2/pi)

运行结果
omega =
18.0263 0 0 0
0 18.0263 0 0
0 0 6.8854 0
0 0 0 6.8854
与给出的结果一致。
例3：采用polyeig函数进行计算
前面我们是将运动方程转换为状态空间之后，通过对状态矩阵进行处理得到固有频率等信息。实际上，Matlab早就考虑了这个问题，对于这样的特征值求解，其提供了一个polyeig函数可以直接进行求解。函数形式很简单，介绍如下：
[转载]Matlab <wbr>状态空间法求解阻尼系统固有频率
那么对于例2中的问题，我们可以直接计算，代码如下：
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;

A2 = M;
A1 = C;
A0 = K;
[X, e] = polyeig(A0, A1, A2)
omega = abs(e / 2/pi)

运行结果：
omega =
18.0263
18.0263
6.8854
6.8854
例4：采用状态空间规范形进行计算
其实对于系统的状态空间描述有很多规范形，比如友形（companion form），模态形（canonical form），而模态形可以直接得到系统的固有频率。
假设原系统的特征值为lambda1, lambda2, alpha + j*beta, alpha - j*beta，则规范形的A矩阵为
[lambda1, 0, 0, 0
0, lambda2, 0, 0
0, 0, alpha, beta
0, 0, -beta, alpha]
matlab提供了canon函数可以直接进行规范形转换。同样的，对例2中的问题进行计算，代码如下。
% 多自由度系统仿真
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;];

% 标准状态方程矩阵

我的更多文章

下载客户端阅读体验更佳

APP专享