ME! 发表于 2013-7-12 16:07

matlab用振型叠加法求不平衡响应

无阻尼系统对不平衡质量的响应,大家帮我看看对不对?f1=m*e我是根据杨橚的机床动力学的公式,质量矩阵已经归一化了,是不是无阻尼的就是这样的啊?

%(3) 正则化特征向量和计算参数
%--------------------------------------------------------------------------
Factor=diag(V1'*M*V1);                           %模态质量
Vnorm=V1*inv(sqrt(diag(Factor)));                  %正则化振型向量(公式5.36即2.2-39)
VnormK=diag(Vnorm'*K*Vnorm);                     %模态刚度
omega2=diag(sqrt(Vnorm'*K*Vnorm)) ;                % 截断固有频率w5.38即2.2-42
% VnormM=diag(Vnorm'*M*Vnorm);                     %正则模态质量   
Fnorm=Vnorm'*ff ;                                     % 模态力矢量

%--------------------------------------------------------------------------
%%(4) 模态坐标的不平衡响应
%--------------------------------------------------------------------------
%eta0=Vnorm'*M*q0; deta0=Vnorm'*M*dq0;            %2.2-53% 初始条件模态坐标的位移和速度               
%eta=zeros(nstep,sdof);                            %%初始化位移2.2-52

phase0=omega0*tt;                                 %w*t omega0为激振力频率

for i=1:sdof                                  % t(i)时刻的响应
gama=omega0/omega(i) ;                         %omega0为激振力的圆频率,omega(i)=w,系统的无阻尼固有频率
X0=1/(1-gama^2);
eta(:,i)=Fnorm(i)*gama^2*X0*(cos(omega0*tt));
end
%--------------------------------------------------------------------------
%   (6) 将模态坐标转化到物理坐标系
%--------------------------------------------------------------------------
eta=eta';            %模态坐标
y=Vnorm*eta;         %物理坐标


%(3) 正则化特征向量和计算参数
%--------------------------------------------------------------------------
Factor=diag(V1'*M*V1);                           %正则化振型
Vnorm=V1*inv(sqrt(diag(Factor)));                  %正则化振型向量(公式5.36即2.2-39)正则化因子alpha=sqrt(diag(Factor)
omega2=diag(sqrt(Vnorm'*K*Vnorm)) ;                % 正则化模态刚度矩阵diag(Vnorm'*K*Vnorm)=w^2 即固有圆频率
VnormM=diag(Vnorm'*M*Vnorm);                           %正则化模态质量矩阵diag(Vnorm'*M*Vnorm)=1
Fnorm=Vnorm'*ff ;                                     % 正则化模态力矢量

%--------------------------------------------------------------------------
%%(4) 模态坐标的脉冲响应
%--------------------------------------------------------------------------
eta0=Vnorm'*M*q0; deta0=Vnorm'*M*dq0;   %2.2-53% 初始条件模态坐标的位移和速度               
eta=zeros(nstep,sdof);                %%初始化位移2.2-52
phase0=omega0*tt;          %w*t omega0为激振力频率

for i=1:sdof                                  % t(i)时刻的响应
gama=omega0/omega(i) ;                         %omega0为激振力的圆频率,omega(i)=w,系统的无阻尼固有频率      
X0=1/(1-gama^2);
eta(:,i)=f1*gama^2*X0*(sin(omega0*tt));
end
%--------------------------------------------------------------------------
%   (6) 将模态坐标转化到物理坐标系
%--------------------------------------------------------------------------
eta=eta';            %模态坐标
y=Vnorm*eta;         %物理坐标
这是y.z方向坐标画的轴心轨迹



简单摇摆 发表于 2013-7-15 20:02

虽然我也在做   但现在还看不懂   纠结   谢谢 啦

dw04116 发表于 2013-7-16 11:40

振型叠加求模态刚度,模态质量是最关键的,你直接省略了。。。
还有你只有一节模态在叠加什么?

ME! 发表于 2013-7-16 16:39

本帖最后由 ME! 于 2013-7-16 16:58 编辑

dw04116 发表于 2013-7-16 11:40 static/image/common/back.gif
振型叠加求模态刚度,模态质量是最关键的,你直接省略了。。。
还有你只有一节模态在叠加什么?
我的模态质量都归一化为1了,那只是某一个自由度方向上的响应模态刚度和模态质量,我求了啊
Factor=diag(V1'*M*V1);                           %模态质量
Vnorm=V1*inv(sqrt(diag(Factor)));                  % 正则化振型向量
VnormK=diag(Vnorm'*K*Vnorm);                     %模态刚度
omega2=diag(sqrt(Vnorm'*K*Vnorm)) ;                % 正则化模态刚度
VnormM=diag(Vnorm'*M*Vnorm);                     %正则化模态质量

下面是叠加之后的程序


请问叠加之后的动态响应图有什么意义,如果要画轴心轨迹不应该是某一自由度的x,y坐标吗?


ME! 发表于 2013-7-16 16:40

本帖最后由 ME! 于 2013-7-16 17:12 编辑

=eig(K,M);                              % 计算特征值和特征向量

Factor=diag(V1'*M*V1);                           %正则化质量矩阵
Vnorm=V1*inv(sqrt(diag(Factor)));                %正则化振型向量(公式5.36即2.2-39)
omega=diag(sqrt(Vnorm'*K*Vnorm)) ;            % 固有频率w^2
VnormM=diag(Vnorm'*M*Vnorm);                      %正则化模态质量

Fnorm=Vnorm'*ff;                                       % 模态力矢量
tt=tt';
=size(K);
=size(tt);
Modamp=Vnorm'*(ac*M+bc*K)*Vnorm;               % 比例阻尼,正则坐标中的阻尼一般是对角阵
zeta=diag((1/2)*Modamp*inv(diag(omega)));                  %   阻尼比 结构的动力特性和响应分析p39(公式2.2-61)
%--------------------------------------------------------------------------
%%(5) 模态坐标的脉冲响应
%--------------------------------------------------------------------------
eta0=Vnorm'*M*q0;
deta0=Vnorm'*M*dq0;                  %2.2-53% 初始条件模态坐标的位移和速度
eta=zeros(nstep,sdof);                %%初始化位移2.2-52
phase0=omega0*tt;                      %w*t omega0为激振力频率
gama=omega0/omega(1);                        %omega0为激振力的圆频率,omega(i)=w,系统的无阻尼固有频率
omegad=omega(1)*sqrt(1-zeta(1)^2);
phase=omegad*tt;                            %wd*t其中(omegad为有阻尼固有频率)
Exx=exp(-zeta(1)*omega(1)*tt);                %中间变量
C1=eta0(1);                                        %初始位移
C2=(deta0(1)+eta0(1)*zeta(1)*omega(1))/omegad;          %中间变量

X0=sqrt((1-gama^2)^2+(2*zeta(1)*gama)^2);
XX=Fnorm(1)/(omega(1)^2*X0);
XP=atan((2*zeta(1)*gama)/(1-gama^2));

D1=(zeta(1)*omega(1)*cos(XP)+omega0*sin(XP))/omegad;
D2=cos(XP);

eta(:,1)=C1*Exx.*cos(phase)+C2*Exx.*sin(phase)...
         -XX*Exx.*(D1*sin(phase)+D2*cos(phase))...
         +XX*cos(phase0-XP);
for i=2:sdof                                  % t(i)时刻的响应
gama=omega0/omega(i);                        %omega0为激振力的圆频率,omega(i)=w,系统的无阻尼固有频率
omegad=omega(i)*sqrt(1-zeta(i)^2);
phase=omegad*tt;                            %wd*t其中(omegad为有阻尼固有频率)
Exx=exp(-zeta(i)*omega(i)*tt);                %中间变量

C1=eta0(i);                                        %初始位移
C2=(deta0(i)+eta0(i)*zeta(i)*omega(i))/omegad;          %中间变量

X0=sqrt((1-gama^2)^2+(2*zeta(i)*gama)^2);
XX=Fnorm(i)/(omega(i)^2*X0);
XP=atan((2*zeta(i)*gama)/(1-gama^2));

D1=(zeta(i)*omega(i)*cos(XP)+omega0*sin(XP))/omegad;
D2=cos(XP);

eta(:,i)= eta(:,i-1)+C1*Exx.*cos(phase)+C2*Exx.*sin(phase)...
         -XX*Exx.*(D1*sin(phase)+D2*cos(phase))...
         +XX*cos(phase0-XP);
                     
end


% 动态响应图
%--------------------------------------------------------------------------
eta=eta';            %模态坐标
y=Vnorm*eta;         %物理坐标

figure(1)
plot(tt,y(jth,:))         %某一自由度的响应
xlabel('时间(seconds)'), ylabel('位移(m)')
title('时间历程响应')

先叠加再返回物理坐标

补充内容 (2013-7-17 16:16):
我觉得这样叠加时错误的

ME! 发表于 2013-7-17 15:45

本帖最后由 ME! 于 2013-7-17 16:15 编辑

y=Vnorm*eta;已经考虑叠加了

ME! 发表于 2013-7-30 19:46

dw04116 发表于 2013-7-16 11:40 static/image/common/back.gif
振型叠加求模态刚度,模态质量是最关键的,你直接省略了。。。
还有你只有一节模态在叠加什么?

能帮我仔细看看吗?真的很需要请教一下!

ME! 发表于 2013-7-30 20:43

我是按照杨橚的《机床动力学》147面的方程画图的,只是变为多自由度

茹果刘儿 发表于 2013-8-1 15:43

过一阵子也要做这个,到时候有问题可以请教你们~~~

1713573225 发表于 2015-12-8 09:00

茹果刘儿 发表于 2013-8-1 15:43
过一阵子也要做这个,到时候有问题可以请教你们~~~

有点看不懂 程序好像不通用吧
页: [1]
查看完整版本: matlab用振型叠加法求不平衡响应