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

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

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

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

x

无阻尼系统对不平衡质量的响应，大家帮我看看对不对？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

登录/注册后可看大图

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

我的模态质量都归一化为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 编辑

[V1,D]=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';
[sdof,n1]=size(K);
[nstep,n2]=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 编辑

QQ截图20130717161348.png

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

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

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

登录/注册后可看大图

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

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

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

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

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

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

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

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

有点看不懂程序好像不通用吧

账号		自动登录	找回密码
密码			我要加入

[转子动力学] matlab用振型叠加法求不平衡响应

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