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方向坐标画的轴心轨迹
虽然我也在做 但现在还看不懂 纠结 谢谢 啦 振型叠加求模态刚度,模态质量是最关键的,你直接省略了。。。
还有你只有一节模态在叠加什么?
本帖最后由 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 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 16:15 编辑
y=Vnorm*eta;已经考虑叠加了
dw04116 发表于 2013-7-16 11:40 static/image/common/back.gif
振型叠加求模态刚度,模态质量是最关键的,你直接省略了。。。
还有你只有一节模态在叠加什么?
能帮我仔细看看吗?真的很需要请教一下! 我是按照杨橚的《机床动力学》147面的方程画图的,只是变为多自由度 过一阵子也要做这个,到时候有问题可以请教你们~~~ 茹果刘儿 发表于 2013-8-1 15:43
过一阵子也要做这个,到时候有问题可以请教你们~~~
有点看不懂 程序好像不通用吧
页:
[1]