当有阻尼时(M——质量矩阵;C—— 阻尼矩阵;K——刚度矩阵)
M*X''+CX'+KX=0
有一个近似算法是现在常用的,就是先求出无阻尼系统的固有频率及相应振型矩阵,用命令eig(K,M,'qz')或者eig(K,M),将系统正则化,求出正则阻尼矩阵,其可能不是对角矩阵,但取其对角元素即可,那么阻尼系统固有频率为wr=w*sqrt(1-c),其中c为正则阻尼矩阵对角线上元素,而系统的振型不变,因此原来地振型矩阵中各列向量仍对应阻尼系统各固有频率
- function [v,w,zeta]=vbr_sf(m,d,k)
- %vbr_sf vbr_sf(m,d,k)
- % [v,w,zeta]=vbr4(m,d,k)
- % function vbr_sf finds the mode shapes and natural frequencies of
- % a linear second order matrix equation.
- % [v,w]=vbr_sf(m,k) finds the mode shapes and natural frequencies
- % for the undamped case.
- if nargin==2
- k=d;
- [v,w]=eig(m\k);
- w=sqrt(w);
- end
- if nargin==3
- if norm(d/m*k-k/m*d) < 1e-8*norm(k/m*d)
- %disp('Damping is proportional, eigenvectors are real.')
- [v,w]=eig(m\k);
- w=sqrt(w);
- zeta=(v'*m*v)\(v'*d*v)/2/w;
- else
- %disp('Damping is non-proportional, eigenvectors are complex.')
- a=[0*k eye(length(k));-m\k -m\d];
- [v,w1]=eig(a);
- w=abs(w1);
- zeta=-real(w1)/w;
- end
- end
- w=diag(w);zeta=diag(zeta);
复制代码
来自:http://blog.sina.com.cn/s/blog_9f81fb09010157c1.html |