复模态分析遇到的问题与经验总结

mxlzhenzhu

一直看书，主要关心比例阻尼系统，但是有的系统阻尼比较大，不合适比例阻尼，遂进一步看复模态的内容。

一直以为自己搞明白了，结果简单列了一下方程，发现一个很有趣的正交性问题还是没有彻底明白，详细内容见附件，恳请论坛大牛一起讨论。





不甚感激。
附件在此： 复模态分析正交性证明.docx (38.82 KB, 下载次数: 16)

2015-3-31 11:29 上传

点击文件名下载附件

mxlzhenzhu

我又问了一个愚蠢的问题，今天发现是因为编程的时候少加了一个.号造成的错误, 造成了正交性问题，结果让人浮想联翩；

MATLAB语法里面，转置和共轭转置是截然不同的概念，今日犯此错误算是刻骨铭心。先把前三种形式的复模态计算传函的MATLAB程序贴出来，感兴趣的可以看看，我的第四种矩阵形式的FRF计算就不再研究了。

1. <p>%% 复模态分析案列
   
2. %% mxl.2015-3-29
   
3. %% 31日，计算结果不理想，暂时不搞了。
   
4. %% 4月8日，共轭转置和转置是有严格区别的;
   
5. clear
   
6. close all
   
7. clc</p><p>
   
8. State_Choice='Second';</p><p>m1=1;m2=10;m3=5.3;m4=0.1;
   
9. c1=10;c2=19;c3=0.8;c4=100;
   
10. k1=1e6;k2=5e6;k3=2.5e5;k4=7e6;
    
11. M=[
    
12.     m1,0,0;
    
13.     0,m2,0;
    
14.     0,0,m3;
    
15. ];
    
16. C=[
    
17.     c1,-c1,0;
    
18.     -c1,c1+c2,-c2;
    
19.     0,-c2,c2+c3;
    
20. ];
    
21. K=[k1,-k1,0;
    
22.     -k1,k1+k2,-k2;
    
23.     0,-k2,k2+k3;
    
24.     ];</p><p>
    
25. % M=[
    
26. %     m1,0,0,0;
    
27. %     0,m2,0,0;
    
28. %     0,0,m3,0;
    
29. %     0,0,0,m4;
    
30. % ];
    
31. % C=[
    
32. %     c1,-c1,0,0;
    
33. %     -c1,c1+c2,-c2,0;
    
34. %     0,-c2,c2+c3,-c3;
    
35. %     0,0,-c3,c4+c3;
    
36. % ];
    
37. % K=[k1,-k1,0,0;
    
38. %     -k1,k1+k2,-k2,0;
    
39. %     0,-k2,k2+k3,-k3;
    
40. %     0,0,-k3,k3+k4;
    
41. %     ];
    
42. N=size(M,1);
    
43. indicei=1;
    
44. indicej=1;
    
45. %   %=====================================================================%
    
46. %   (1)State transform equation.
    
47. %   %=====================================================================%
    
48. if strcmp(State_Choice,'First')
    
49.     A=[
    
50.         C,K;
    
51.         K,zeros(size(M));
    
52.     ];</p><p>    B=[
    
53.         M,zeros(size(M));
    
54.         zeros(size(M)),-K;
    
55.         ];% I invented
    
56.     indiceii=indicei;
    
57.     indicejj=indicej+N;</p><p>elseif strcmp(State_Choice,'Second')
    
58.     A=[
    
59.         -M,zeros(size(M));
    
60.         zeros(size(M)),K;
    
61.         ];
    
62.    
    
63.     B=[
    
64.         zeros(size(M)),M;
    
65.         M,C;
    
66.     ];% Ji wen mei
    
67.     indiceii=indicei+N;
    
68.     indicejj=indicej+N;
    
69. elseif strcmp(State_Choice,'Third')
    
70.     A=[
    
71.         K,zeros(size(M));
    
72.         zeros(size(M)),-M;
    
73.         ];
    
74.    
    
75.     B=[
    
76.         C,M;
    
77.         M,zeros(size(M))];% D.J.Ewins
    
78.     indiceii=indicei;
    
79.     indicejj=indicej;   
    
80. else
    
81.     B=[
    
82.         M,C;
    
83.         zeros(size(M)),M;
    
84.     ];
    
85.     A=[
    
86.         zeros(size(M)),K;
    
87.         -M,zeros(size(M));
    
88.         ];% my word
    
89.     indiceii=indicei;
    
90.     indicejj=indicej+N;
    
91. end
    
92. %   %=====================================================================%
    
93. %   (2)Eigen Equ.
    
94. %   %=====================================================================%
    
95. [V,D]=eig(A,-B);
    
96. nameda=diag(D);
    
97. % return
    
98. %   %=====================================================================%
    
99. %   (3)FRF, modal superpostion
    
100. %   %=====================================================================%</p><p>frequency=[2:2:800]';
     
101. H=cell(size(frequency));
     
102. % 基于复振型和复特征值计算FRF
     
103. for loopi=1:numel(frequency)
     
104.     omega=2*pi*frequency(loopi);
     
105.     for loop=1:numel(nameda)
     
106.         ar=V(:,loop).'*B*V(:,loop);
     
107.         if loop>1
     
108.             H{loopi,1}=H{loopi,1}+V(:,loop)*V(:,loop).'/ar/(1j*omega-nameda(loop));
     
109.         else
     
110.             H{loopi,1}=V(:,loop)*V(:,loop).'/ar/(1j*omega-nameda(loop));
     
111.         end
     
112.     end
     
113. end</p><p> </p><p> </p><p> </p><p>FRF=zeros(size(frequency));
     
114. for loop=1:numel(frequency)
     
115.     FRF(loop,1)=H{loop,1}(indiceii,indicejj);
     
116. end</p><p>figure
     
117. subplot(2,1,1)
     
118. plot(frequency,abs(FRF),'r')
     
119. subplot(2,1,2)
     
120. plot(frequency,angle(FRF),'k')
     
121. %   %=====================================================================%
     
122. %   (4)FRF, direct
     
123. %   %=====================================================================%
     
124. omega=2*pi*frequency;
     
125. H=cell(size(frequency));
     
126. for loop=1:numel(frequency)
     
127.     H{loop,1}=inv(K-omega(loop)^2*M+1j*omega(loop)*C);
     
128. end</p><p> </p><p>FRF=zeros(size(frequency));
     
129. for loop=1:numel(frequency)
     
130.     FRF(loop,1)=H{loop,1}(indicei,indicej);
     
131. end</p><p>figure
     
132. subplot(2,1,1)
     
133. plot(frequency,abs(FRF),'r')
     
134. subplot(2,1,2)
     
135. plot(frequency,angle(FRF),'k')
     
136. </p>

复制代码

mxlzhenzhu

Two wrongs make a right
{:{03}:}{:{03}:}

zswseu

感觉像ERA中的一段。
很多理论都很深奥的。