发布功能喝代码高亮

Posted 2022-05-19 linyunfeng2019

ss模型复模态的物理意义及adams复模态振型求解
2017年03月28日 19:48:37 lijil168 阅读数 453
adams可以线性化输出SS状态空间ABCD矩阵，根据A可以进行系统复模态计算。

A可看作是角频率张量矩阵，由模态矩阵进行坐标变换解耦后可以得到所有的模态参数，特征向量（模态系数）为复数不仅影响幅值，还影响相位角，而特征值为复数，其虚部为阻尼频率，而实部会使振动随时间衰减。 

 

GX˙+HX=0
ФTGФ  Ф-1X˙+ФTHФ  Ф-1X=0
将X˙转换到主模态基下，G转换到主模态对偶基下，在对偶基下完成系统解耦。
[en]T ФTGФ  Ф-1X˙+ [en]T ФTHФ Ф-1X=0
 
 
X˙=-G\HX=AX
A可认为是角频率张量矩阵。
Ф-1X˙=Ф-1 (ФT )-1 ФTAФ Ф-1 X
将位移X和速度X˙由逆模态振型矩阵Ф-1变换到主坐标基中，角频率阵经ФTAФ变换到主基的对偶基后，乘以位移矩阵，得到对偶基下的速度矩阵，再经Ф-1 (ФT )-1变换到主基下，在主基下完成系统解耦。
以张量描述为：
[ep]TФ-1X˙= [ep]TФ-1 (ФT)-1 ФTAФ Ф-1 X
 
设位移矢量、速度矢量、角频率张量在原基下：Er :[e1 e2]T为单位基矢量。
位移矢量为：[e1 e2][x1 x2]T,
速度矢量为：[e1 e2][x1˙ x2˙]T,
角频率张量为：[e1 e2][A11, A12;A21,A22] [e1 e2] T。
 
主基为Ep: [ep1 ep2]T为单位基矢量。
主基的对偶基为En: [en1 en2]T为单位基矢量。
 
Er到Ep的过渡矩阵为：Arp=Ф=[e1 e2]T [ep1 ep2]
Ep到Er的过渡矩阵为：Apr=Ф-1=[ep1 ep2]T [e1 e2]
En到Er的过渡矩阵为：Anr=ФT=[en1 en2]T [e1 e2]
En到Ep的过渡矩阵为：Anp=ФTФ=[en1 en2]T [ep1 ep2]
Ep到En的过渡矩阵为：Apn=(ФTФ)-1=[ep1 ep2]T [en1 en2]
 
带入上述SS状态方程可得：
[ep1 ep2] [ep1] [e1 e2]  [x1˙]   =
[ep2]         [ x2˙]
   
[ep1 ep2]   [ep1] [en1 en2]  [en1 ] [e1 e2]  [A11, A12]   [e1]  [ep1 ep2]  [ep1] [e1 e2] [x1]
[ep2]         [en2]         [A21,A22]  [e2]         [ep2]         [x2]
 
从而可以清楚地看到坐标变换（解耦）的过程。

A=[-1.00377172158969996E+00 -3.91715793777936216E+00  0.00000000000000000E+00 -1.96627966300476429E-06
 1.00000000000000000E+00 0.00000000000000000E+00  0.00000000000000000E+00  0.00000000000000000E+00
 0.00000000000000000E+00 5.65305403113869734E-06 -1.00377172158969996E+00 -1.12618290711156668E+01
 0.00000000000000000E+00 0.00000000000000000E+00 1.00000000000000000E+00  0.00000000000000000E+00
];
%A=load(‘f:\temp\tta.txt‘);
 f_M_K=A([1,3],[2,4]);
 M_K=-f_M_K;
 [vec,ome]=eig(M_K);
 omeg=diag(sqrt(ome));
 omeg_hz=omeg/2/pi;
 f_M_C=A([1,3],[1,3]);
 M_C=-f_M_C;
 ssA=[zeros(2),eye(2);-M_K,-M_C];
 [Vec,Ome]=eig(ssA);
 Omeg=imag(diag(Ome([1,3],[1,3])));
 alpha=real(diag(Ome([1,3],[1,3])));
 Omeg_hz=Omeg/2/pi;  %%自然振荡频率
 alpha_hz=alpha/2/pi;
 Vec_single=Vec(:,[1,3]);
 Vec_abs=abs(Vec_single);
 Vec_ang=angle(Vec_single)*180/pi;
 
 mode_vec=vec;
 mode_vec(:,1)=mode_vec(:,1)./mode_vec(1,1);
mode_vec(:,2)=mode_vec(:,2)./mode_vec(1,2);
subplot(3,3,1)
plot([1;2],mode_vec(:,1))
title(strcat(‘w1=‘,num2str(omeg_hz(1))));
subplot(3,3,2)
plot([1;2],mode_vec(:,2))
title(strcat(‘w2=‘,num2str(omeg_hz(2))));
 
 
 mode_abs=Vec_abs(1:2,:);
 mode_ang=Vec_ang(1:2,:);
 mode_abs(:,1)=mode_abs(:,1)./mode_abs(1,1);
mode_abs(:,2)=mode_abs(:,2)./mode_abs(1,2);
subplot(3,3,3)
plot([1;2],mode_abs(:,1))
title(strcat(‘alpha1=‘,num2str(alpha_hz(1)),‘beta1=‘,num2str(Omeg_hz(1))));
subplot(3,3,4)
plot([1;2],mode_abs(:,2))
title(strcat(‘alpha2=‘,num2str(alpha_hz(2)),‘beta2=‘,num2str(Omeg_hz(2))));

subplot(3,3,5)
plot([1;2],mode_ang(:,1))
title(strcat(‘ang_mode1‘,‘beta1=‘,num2str(Omeg_hz(1))));

subplot(3,3,6)
plot([1;2],mode_ang(:,2))
title(strcat(‘ang_mode2‘,‘beta2=‘,num2str(Omeg_hz(2))));

subplot(3,3,7)
polar(mode_ang(:,1)*pi/180,mode_abs(:,1));
title(strcat(‘ang_mode1‘,‘beta1=‘,num2str(Omeg_hz(1))));
subplot(3,3,8)
polar(mode_ang(:,2)*pi/180,mode_abs(:,2));
title(strcat(‘ang_mode2‘,‘beta2=‘,num2str(Omeg_hz(2))));

  目前好像不支持matlab代码高亮

for j = 1:NumStates
    if imag(l(j))~=0
       %scale the complex eigenvectors so that the maximum element is 1+j0
       [maxu,mu_idx] = max(abs(u(:,j)));
       u(:,j) = u(:,j)/u(mu_idx,j);
    end
end