声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 6813|回复: 22

[编程技巧] 通过状态方程求固有频率、振型遇到的问题

  [复制链接]
发表于 2012-4-17 10:28 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 不会不怕 于 2012-4-17 10:28 编辑

M=diag([Mb,Jb,Mf1,Jf1,Mf2,Jf2]);
K=[2*Ks,0,-Ks,0,-Ks,0;...
   0,2*Ks*(lb^2),Ks*lb,0,-Ks*lb,0;...
   -Ks,Ks*lb,Ks+2*Kp,0,0,0;...
   0,0,0,2*Kp*(lf^2),0,0;...
   -Ks,-Ks*lb,0,0,Ks+2*Kp,0;...
   0,0,0,0,0,2*Kp*(lf^2)];
C=[2*Cs,0,-Cs,0,-Cs,0;...
   0,2*Cs*(lb^2),Cs*lb,0,-Cs*lb,0;...
   -Cs,Cs*lb,Cs+2*Cp,0,0,0;...
   0,0,0,2*Cp*(lf^2),0,0;...
   -Cs,-Cs*lb,0,0,Cs+2*Cp,0;...
   0,0,0,0,0,2*Cp*(lf^2)];
B=eye(size(M));
F=zeros(size(M));
E=inv(M);
G=-E*K;
H=-E*C;
A=[F,B;G,H];
% 3. 求矩阵的特征值与特征向量
[V,D]=eig(A);  % 矩阵A的特征值(D)与特征向量(V)
% 4. 计算相应参数abs(imag(D(i,i)))
for i=1:size(A)
    omega_n(i)=sqrt(real(D(i,i))^2+imag(D(i,i))^2*i);
    omega_d(i)=abs(imag(D(i,i)));
    sigama(i)=sqrt(1-(omega_d(i)/(omega_n(i))^2))
end
这个程序求出的固有频率不正确,不知道哪里出现了问题,请各位指导,谢谢
回复
分享到:

使用道具 举报

发表于 2012-4-17 10:57 | 显示全部楼层
回复 1 # 不会不怕 的帖子

A=[C M;M 0];B=[K 0;0 -M];
[V,D]=eig(-B,A)
看看这些能不能解决你的问题。

评分

1

查看全部评分

 楼主| 发表于 2012-4-18 09:51 | 显示全部楼层
回复 2 # yyxt007 的帖子

这样的结果还是不正确,谢谢,其实我觉得是不是状态方程后通过特征向量求固有频率的公式就不对了呢,一直没找到方法
发表于 2012-4-18 09:59 | 显示全部楼层
回复 3 # 不会不怕 的帖子

你所说的不正确,指的是什么?
 楼主| 发表于 2012-4-18 10:24 | 显示全部楼层
回复 4 # yyxt007 的帖子

23.5254        23.5254        4.0759        4.0759        20.9088        20.9088        34.9359        34.9359        6.6887        6.6887        34.9359        34.9359
这是matlab运行出来的有阻尼固有频率,与结构实际频率不符,实际频率在10hz左右
发表于 2012-4-18 10:33 | 显示全部楼层
回复 5 # 不会不怕 的帖子

最好把频率重新排列;
计算结果与实际频率不符,是程序的问题还是算法不对?
实际结构的频率,是测量得到的?
对于一般结构,计算频率是不需考虑阻尼的。
我用过我提到的代码,计算结果没有问题。

评分

1

查看全部评分

 楼主| 发表于 2012-4-18 15:00 | 显示全部楼层
回复 6 # yyxt007 的帖子

是简化的车辆振动模型,将矩阵方程式,改写成状态方程之后计算的,与车辆实际频率不符。可能您的代码我运算错了
发表于 2012-11-5 00:08 | 显示全部楼层
如果不用状态方程方法,固有频率应该是eig(inv(M)*K),我说的这样算出来的特征值肯定是系统的固有频率的。可以再振动力学上,多自由度固有频率中提到。至于状态方程方法算出来的特征值可能不是固有频率

点评

状态方程方法算出来的特征值是固有频率  发表于 2012-11-6 22:38
发表于 2012-11-5 00:09 | 显示全部楼层
我已在找,怎么用状态方程求固有频率,求特征值是找到了,不知道是不是等价的
发表于 2012-11-7 12:48 | 显示全部楼层
问题在与两种方法算出来系统的特征值明显不一样,状态空间是eig([A])。(其中B=eye(size(M));
F=zeros(size(M));
E=inv(M);
G=-E*K;
H=-E*C;
A=[F,B;G,H];),
而振动 多自由度分析是eig(inv(M)*K)。当然,一般这个不考虑阻尼,
两钟方法结果以那个为准??
发表于 2012-11-7 23:31 | 显示全部楼层
太久没玩这有些生疏了
MX"+KX=0
x1=X, x2=X' => x1'=X'=x2 & x2'=X"=-inv(M)*K*x1 => aa=[0,I;-inv(M)*K 0]
  1. clc; clear;
  2. M=[2 0;0 4]; K=[8 0;0 64]; eig(inv(M)*K)
  3. aa=[zeros(size(M)) eye(size(M));-inv(M)*K zeros(size(M))]; eig(aa)
复制代码
别忘记前面是在求w^2
回复 支持 0 反对 1

使用道具 举报

发表于 2012-11-30 09:33 | 显示全部楼层
这类问题建议从模型本身开始入手,首先要确定模型本身是否正确,否则一切都是空谈
发表于 2012-12-9 19:31 | 显示全部楼层
请问非对角阵系数的方程怎么求响应,因为列微分方程时就不止一个二阶导数的未知数了,
还有就是怎么用状态方程解
发表于 2012-12-11 10:21 | 显示全部楼层

1. 对你的表述不是很理解,最好能给出数学模型来
2. 你的第二个问题就是将常见的高阶微分方程组,转换到状态空间中,变成一阶状态空间方程,然后求解,具体的变换方法论坛很多帖子都有,可以搜索一下
发表于 2012-12-11 10:33 | 显示全部楼层
mx''+kx=f(t),m为非对角阵,可以列出n个微分方程,例如m(1,1)x''+m(1,4)x''+k(1,1)x+k(1,1)x=f(t)
我看见matlab里面可以用ode23解微分方程,但是我的微分方程,我不知道怎么代换才可以变换成matlab能解的形式
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 01:59 , Processed in 0.084288 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表