如何通过状态方程求固有频率

ChaChing

昨晚看/回帖, 学习并复习了一下多年前的一些回覆!

求一个求四自由度振动系统固有频率的MATLAB程序 http://forum.vibunion.com/thread-70695-1-1.html
通过状态方程求固有频率、振型遇到的问题 http://forum.vibunion.com/thread-110946-1-1.html
多自由度非比例阻尼的线性方程解法！！ http://forum.vibunion.com/thread-44115-1-1.html

由上头这些帖不难发现有因对比不正确, 就怀疑通过状态方程的特徵矩阵算出来的特征值可能不是固有频率!? 请别再质疑了, 一定是没注意到某些细节! google下网上资料或者查看下Ogata的书(Modern Control Engineering)吧!
或者到处求助, 但可能没经验, 忽略掉给出细节没能描述明确! 没头没尾仅是简单叙述, 至少在个人有限水平下, 仅能确认somthing is wrong, 但不知从那帮忙debug起!?

想想好像只能自己整理一下, 希望遇到问题的人可以自己找出差异, 更正自己原先的错误!
应该不是很难才是, 但个人真有些力不从心且也很懒, 就儘可能交代一下了! 希望能说明白

首先稍微推导下, 一些相对应的式子及与matlab如何对应, 由於个人水平有限, 仅为了后面自己撰写程序方便而已, 应该不是很清楚严谨! 凑合参考看吧, 如果想要详细了解(如高阶或其它型式...), 就请看下一些相关的书吧!

从简易起, 无阻尼系统  Mx"+Kx=0

1. x=q*exp(i*w) => (-w^2*M+K)*q=0 => K*q=w^2*M*q
    matlab eig函数 A*x=lambda*B*x => A=K; B=M; lambda=w^2;
2. x=q*exp(i*w) => (-w^2*M+K)*q=0 => inv(M)*K*q=w^2*q
    matlab eig函数 A*x=lambda*x => A=inv(M)*K; lambda=w^2; 或 A=M\K; lambda=w^2;
3. x1=x, x2=x' => x1'=x'=x2 & x2'=x"=-inv(M)*K*x1
    => 特徵矩阵 A=[0,I;-inv(M)*K 0] 或  A=[0,I;-M\K 0]
4. eqn 4.6.11 of http://forum.vibunion.com/thread-44115-1-1.html
   => A=[0 M;M 0]; B=[K 0;0 -M]; => lambda*A*q=-B*q (eqn 4.6.14)

注意下,方式1/2对应一般方式, 方式3/4对应状态方程方式, 其中inv函数可以用倒除取代(算法不同,速度比较快)
另提醒下状态方程是有许多种形式的, 其特徵矩阵的特徵值会形成共軛复根, 其虚部即为w(=2*pi*f), 参见Ogata的书

有阻尼系统 Mx"+Cx'+Kx=0

1. x1=x, x2=x' => x1'=x'=x2 & x2'=x"=-inv(M)*K*x1-inv(M)*C*x2
   => A=[0,I;-inv(M)*K -inv(M)*C] 或  A=[0,I;-M\K -M\C]
2. eq4.6.11 of http://forum.vibunion.com/thread-44115-1-1.html
   => A=[C M;M 0]; B=[K 0;0 -M]; => lambda*A*q=-B*q (eqn 4.6.14)

ChaChing

为试下各方式的结果差异, 先造个简易型的M/K矩阵来试吧!
两个自由度故意没耦合(仅对角项有值), 如此可以预知其两个特征频率为w1=2及w2=4(注意下w为角频率)
有阻尼系统中的C矩阵值, 也特意取成两个自由度的阻尼比皆为0.05

无阻尼系统  Mx"+Kx=0

1. clc; clear;
   
2. M=[2 0;0 4]; K=[8 0;0 64]; Z0=zeros(size(M)); I1=eye(size(M));
   
3. A=K; B=M; lambda=eig(A,B); w1=sqrt(lambda);
   
4. A=inv(M)*K; lambda=eig(A); w2a=sqrt(lambda);
   
5. A=M\K; lambda=eig(A); w2b=sqrt(lambda);
   
6. A=[Z0 I1;-inv(M)*K Z0]; lambda=eig(A); w3a=imag(lambda);
   
7. A=[Z0 I1;-M\K Z0]; lambda=eig(A); w3b=imag(lambda);
   
8. A=[Z0 M;M Z0]; B=[K Z0;Z0 -M]; lambda=eig(-B,A); w4=imag(lambda);
   
9. [w1 w2a w2b], [w3a w3b w4]
   

复制代码

[w1 w2a w2b] =
    2.0000    2.0000    2.0000
    4.0000    4.0000    4.0000

[w3a w3b w4] =
    2.0000    2.0000    2.0000
   -2.0000   -2.0000   -2.0000
    4.0000    4.0000    4.0000
   -4.0000   -4.0000   -4.0000

有阻尼系统 Mx"+Cx'+Kx=0

1. clc; clear;
   
2. M=[2 0;0 4]; C=[0.4 0;0 1.6]; K=[8 0;0 64]; dapR=0.05; Z0=zeros(size(M)); I1=eye(size(M));
   
3. A=K; B=M; lambda=eig(A,B); w1=sqrt(lambda)*sqrt(1-dapR*dapR);
   
4. A=[Z0 I1;-inv(M)*K -inv(M)*C]; lambda=eig(A); w3a=imag(lambda);
   
5. A=[Z0 I1;-M\K -M\C]; lambda=eig(A); w3b=imag(lambda);
   
6. A=[C M;M Z0]; B=[K Z0;Z0 -M]; lambda=eig(-B,A); w4=imag(lambda);
   
7. w1, [w3a w3b w4]

复制代码

w1 =
    1.9975
    3.9950

[w3a w3b w4] =
    1.9975    1.9975    1.9975
   -1.9975   -1.9975   -1.9975
    3.9950    3.9950    3.9950
   -3.9950   -3.9950   -3.9950

本来还想试一下, "求一个求四自由度振动系统固有频率的MATLAB程序" http://forum.vibunion.com/thread-70695-1-1.html中比较实际的例子, 但刚发现该楼主没给齐相关资料也懒了, 就不试了
有兴趣的可以自行试下有非对角项的例子

牛小贱

先收藏再说，貌似很有用！

linda8866

这么好的帖子必须顶。谢谢楼主，感谢分享。

mayaview

ChaChing 发表于 2013-11-26 16:11
为试下各方式的结果差异, 先造个简易型的M/K矩阵来试吧!
两个自由度故意没耦合(仅对角项有值), 如此可以预 ...

不明觉厉啊！状态方程在哪里？

随心就动

很专业啊，学习了

springen

学习了！