线性加速度法和模态叠加计算振动共振点不重合
上面的图片是通过newmark法计算的单频率激励力下的梁的某点的振动响应时域结果通过matlab内置的fft变换后得到的频域图。下面的图片是通过分析特征值得到固有特性，然后用模态叠加法计算得到响应图。
看共振频率，一阶模态的响应最大，但是两种分析结果的共振点略有差别，不知道是什么原因造成的，是理论本身，还是我的计算过程。
如果考虑阻尼影响，共振频率一般会比固有频率小一些。newmark法计算会有数值阻尼影响。 回复 2 # yyxt007 的帖子
你说的没错，
但是我的模态叠加的原理中也是用模态刚度，模态阻尼，模态质量这些计算的，也就是说也是考了峰值频率不同于固有频率的问题，我这里就相当于比较峰值频率，依然有点差别。。单自由度系统则没有这个差别…… 回复 1 # dw04116 的帖子
newmark方法是从短时间时域上的分析。分析结果与截取时间段有关。
而振型叠加法，如果直接由振型矢量上得到模态，则其具有统计特征，是一个宏观的结果。
请问楼主我用振型叠加法求出了时间响应图，然后进行fft变换，得到的频率响应图的共振频率特别小时怎么回事？激励频率为1500HZ ME! 发表于 2013-6-26 15:40 static/image/common/back.gif
请问楼主我用振型叠加法求出了时间响应图，然后进行fft变换，得到的频率响应图的共振频率特别小时怎么回事？ ...
请教，振型叠加法如何计算时域图 本帖最后由 ME! 于 2013-7-1 15:51 编辑
dw04116 发表于 2013-7-1 11:17 static/image/common/back.gif
请教，振型叠加法如何计算时域图
有公式啊，参考徐斌《matlab有限元结构动力学分析及工程应用》那你用振型叠加法计算的频率响应图是怎么计算出来的，不是利用时域响应图FFT变换的来的吗？
dw04116 发表于 2013-7-1 11:17 static/image/common/back.gif
请教，振型叠加法如何计算时域图
请问你用newmark法求响应的时候得出的共振频率正好是你计算的固有频率吗？我怎么计算的时候得到的幅值最大的点不是对应固有频率的点啊，我看徐荣桥的《结构分析的有限元法与matlab程序设计》里面对初始条件的响应对应的正好是模态分析计算得到的固有频率点 固有频率没问题，但是时域计算时，固有频率处的共振有一个能量累加的过程，所以就不是很吻合了 ME! 发表于 2013-7-1 15:32 static/image/common/back.gif
有公式啊，参考徐斌《matlab有限元结构动力学分析及工程应用》那你用振型叠加法计算的频率响应图是怎么计 ...
不是哎。。不能一个fft就叫振型叠加，它只是把时域信号变成频域而已。我的振型叠加是通过求特征值然后是特征矩阵，再模态刚度、模态质量，然后叠加计算的 本帖最后由 ME! 于 2013-7-4 15:32 编辑
dw04116 发表于 2013-7-4 12:35 static/image/common/back.gif
不是哎。。不能一个fft就叫振型叠加，它只是把时域信号变成频域而已。我的振型叠加是通过求特征值然后是特 ...
你是说我计算的是节点某自由度的时域响应，不是系统的，系统的是要所有自由度的时域响应叠加然后fft求出频率响应对吧，我昨天看了书意识到这个问题了
徐斌的《matlab 有限元结构动力学与工程应用》这本书里面只是某个自由度的时域响应，要求系统的响应得所有自由度的叠加
可是我用newmark计算时不应该会出现频率特别小不等于固有频率的情况啊，不知道为什么计算出来的共振频率也特别小
补充内容 (2013-7-5 14:56):
可是我叠加之后频率还是很小，不等于固有频率
补充内容 (2013-7-8 21:06):
能不能告诉我下你的联系方式啊，能即时联系你的，我想把这个问题弄清楚！老是碰不见你啊！ 本帖最后由 ME! 于 2013-7-8 21:48 编辑
下面是我的程序，能帮我看看吗，共振处的频率不等于omega1，我用振型叠加法也是得不到共振频率，不知道是哪里错了%-------------------------------------------------------------------------%
%--------------------------------------------------------------------------
%newmark
%(0) 输入控制数据
%--------------------------------------------------------------------------
clear; clc;
jk=3;
c=; %约束的节点两点支撑，4为前轴承为节点，13为后轴承位节点
jth=14; % 第jth自由度方向的响应
=size(c);
No_el=16; % 单元数
No_nel=2; % 每个单元的节点数
No_dof=3; % 每个节点自由度数
No_node=No_el+1; % 系统的节点数
Sys_dof=No_node*No_dof; % 系统的自由度数
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
ac=0.00002;bc=0.00000000008; % 比例阻尼参数
al=0; % 空间梁单元整体坐标系和局部坐标系之间的角度
omega0=1600; % 施加的激振力的频率
%(2) 施加节点载荷
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%(3) 计算单元矩阵和组装
L=;%长度分段
d1= ;%主轴分段外径
d2= ;%主轴分段内径
dm1=;%主轴分段附加零件外径
dm2=;%主轴分段附加零件内径
K=sysK(L,d1,d2); %组装刚度矩阵
M=sysM(L,d1,d2,dm1,dm2); %组装质量矩阵
yueshu = [ c(1), 1, 0; c(2),1, 0 ];%约束第c(1)、c(1)节点的轴向位移为0
Kr1=100e6;%前轴承的径向刚度
Kr2=100e6;%后轴承的径向刚度
Kr3=100e6;%中间轴承的径向刚度
=combination_surface_stiffness(Kr1,Kr2,Kr3); %计算轴承结合面刚度
K(c(1)*3-1,c(1)*3-1)=K(c(1)*3-1,c(1)*3-1)+KK1;%将前轴承径向结合部刚度值并入整体刚度矩阵
K(c(2)*3-1,c(2)*3-1)=K(c(2)*3-1,c(2)*3-1)+KK2;%将后轴承径向结合部刚度值并入整体刚度矩阵
= size( yueshu ) ; %取出约束节点矩阵的行数
w2max = max( diag(K)./diag(M) ) ; %修改整体质量矩阵的中间变量
%--------------------------------------------------------------------------
%（4）处理第一类约束条件，修改刚度矩阵和质量矩阵。（采用划行划列法）
for ibc=1:1:bc1_number
cc = yueshu(ibc, 1 ) ;
d = yueshu(ibc, 2 ) ;
m = (cc-1)*3 + d ;
K(:,m) = zeros(No_node*3, 1 ) ;
K(m,:) = zeros( 1, No_node*3 ) ;
K(m,m) = 1;
M(:,m) = zeros( No_node*3, 1 ) ;
M(m,:) = zeros( 1, No_node*3 ) ;
M(m,m) = K(m,m)/w2max/1e5 ;
end
= eigs(K, M,jk,'SM') ; %SM,计算前3阶最小的特征值
%（5）修改特征向量中受约束的自由度
for ibc=1:1:bc1_number
cc = yueshu(ibc, 1 ) ;
d = yueshu(ibc, 2 ) ;
m = (cc-1)*3 + d ;
gEigVector(m,:) = yueshu(ibc,3) ;
end
omega1=sort(sqrt(diag(gEigValue))/(2*pi)) % %约束状态下的固有频率
%--------------------------------------------------------------------------
%(6) 计算瞬态响应
%--------------------------------------------------------------------------
cc1=ac*M+bc*K; % 比例阻尼矩阵
gDeltaT = 0.01 ; %步长
gTimeEnd = (2^14)*gDeltaT; % 计算时间为载荷通过所需时间的两倍
timestep = floor(gTimeEnd/gDeltaT) ; %载荷步
tt=0:0.01:gTimeEnd;
ff=zeros(Sys_dof,1); % 初始化激振力矢量
P=; % 施加激振力的幅值
ff(No_dof*(P(2)-1)+P(3))=P(1); %将力施加在某个自由度上
u=cos(omega0*tt); %简谐激振力
ff=ff*u;
% 定义位移，速度和加速度
gDisp = zeros( (No_el+1)*3, timestep ) ;
gVelo = zeros( (No_el+1)*3,timestep ) ;
gAcce = zeros( (No_el+1)*3, timestep ) ;
% 初始条件
gDisp(:,1) = zeros( (No_el+1)*3, 1 ) ;
gVelo(:,1) = ones( (No_el+1)*3, 1 ) ;
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% 动态响应图
%--------------------------------------------------------------------------
% step3.1 初始计算
gama = 0.5 ;
beta = 0.25 ;
C = zeros( size( K ) ) ;
= size( K ) ;
alpha0 = 1/beta/gDeltaT^2 ;
alpha1 = gama/beta/gDeltaT ;
alpha2 = 1/beta/gDeltaT ;
alpha3 = 1/2/beta - 1 ;
alpha4 = gama/beta - 1 ;
alpha5 = gDeltaT/2*(gama/beta-2) ;
alpha6 = gDeltaT*(1-gama) ;
alpha7 = gama*gDeltaT ;
K1 = K + alpha0*M + alpha1*C;
timestep = floor(gTimeEnd/gDeltaT) ;
% step3.2 对K1进行边界条件处理
= size( yueshu ) ;
K1im = zeros(N,bc1_number) ;
for ibc=1:1:bc1_number
n = yueshu(ibc, 1 ) ;
d = yueshu(ibc, 2 ) ;
m = (n-1)*3 + d ;
K1im(:,ibc) = K1(:,m) ;
K1(:,m) = zeros( No_node*3, 1 ) ;
K1(m,:) = zeros( 1, No_node*3 ) ;
K1(m,m) = 1.0 ;
end
= lu(K1) ; % 进行三角分解，节省后面的求解时间
% step3.3 计算初始加速度
Acce(:,1) = M\(ff(:,1)-K*gDisp(:,1)-C*gVelo(:,1)) ;
% step3.4 对每一个时间步计算
for i=2:1:timestep
fprintf( '当前时间步：%d\n', i ) ;
f1 =ff(:,i)+M*(alpha0*gDisp(:,i-1)+alpha2*gVelo(:,i-1)+alpha3*gAcce(:,i-1)) ...
+ C*(alpha1*gDisp(:,i-1)+alpha4*gVelo(:,i-1)+alpha5*gAcce(:,i-1)) ;
% 对f1进行边界条件处理
= size( yueshu ) ;
for ibc=1:1:bc1_number
n = yueshu(ibc, 1 ) ;
d = yueshu(ibc, 2 ) ;
m = (n-1)*3 + d ;
f1 = f1 - yueshu(ibc,3) * K1im(:,ibc) ;
f1(m) = yueshu(ibc,3) ;
end
y = KL\f1 ;
gDisp(:,i) = KU\y ;
gAcce(:,i) = alpha0*(gDisp(:,i)-gDisp(:,i-1)) - alpha2*gVelo(:,i-1) - alpha3*gAcce(:,i-1) ;
gVelo(:,i) = gVelo(:,i-1) + alpha6*gAcce(:,i-1) + alpha7*gAcce(:,i) ;
end
%--------------------------------------------------------------------------
%时域的快速傅里叶变换 y(n,:)
%--------------------------------------------------------------------------
% 绘制时程曲线
% = size(gNode) ;
t = 0:gDeltaT:gTimeEnd-gDeltaT ;
d = gDisp((floor(No_node/4)*3)+2,:) ;
v = gVelo((floor(No_node/4)*3)+2,:) ;
a = gAcce((floor(No_node/4)*3)+2,:) ;
% tt = 0:gDeltaT/10:gTimeEnd-gDeltaT ;
fd = fft( d ) ;
df = 1/gTimeEnd ;%FFT变换的频率分辨率
f = (0:length(d)-1)*df ;
figure
subplot(211);plot(t,d)
subplot(212);plot(f,abs(fd)) ;
title( 'L/4处挠度的频谱图' ) ;
set(gcf,'color',);
%set(0,'defaultaxesfontsize',8)
xlabel( '频率(Hz)') ;
ylabel( '幅度' ) ;
附件是我的频谱图和响应图 本帖最后由 ME! 于 2013-7-9 20:27 编辑
振型叠加法dt=0.0001; % 时间步长
ti=0; % 初始时间
tf=1; % 结束时间
nt=fix((tf-ti)/dt); % 时间载荷步
tt=ti:dt:ti+nt*dt; % 时间矢量
=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
Fnorm=Vnorm'*ff; % 模态力矢量
tt=tt';
=size(K);
=size(tt);
Modamp=Vnorm'*(ac*M+bc*K)*Vnorm; % 比例阻尼，正则坐标中的阻尼一般是对角阵
zeta=diag((1/2)*Modamp*inv(diag(omega2))); % 阻尼比 结构的动力特性和响应分析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'; %模态坐标
yt=Vnorm*eta; %物理坐标
figure(1)
plot(tt,yt)
xlabel('时间(seconds)'), ylabel('位移(m)')
title('时间历程响应')
tt=tt';yt=yt';
ntime=max(size(tt)); %数据点数
dt=(tt(1,ntime)-tt(1,1))/(ntime-1);%步长
N=fix(log10(ntime)/log10(2));%(2) 截断数据点 y
yfft=fft(yt(1:2^N,:));%(3) 快速傅里叶变换（取绝对值）
yfft=abs(yfft(1:2^N/2,:))*dt;%幅值的绝对值
freq0=0;
freqf= (1/dt)/2; %最大频率值
df=freqf/(2^N/2); %频率间隔
freq=0:df:freqf-df; %频率值
figure (2)
plot(freq,yfft) ;
title( '频谱图' ) ;
xlabel( '频率(Hz)') ;
ylabel( '幅度' ) ;得出来的频谱图也不对！ ME! 发表于 2013-7-9 20:11 static/image/common/back.gif
振型叠加法得出来的频谱图也不对！
很难一一去看。
我感觉的你搞的太复杂了，就算是想用有限元法也没有必要去计算单元矩阵，直接套用现有的单元矩阵多好啊。然后模态叠加法，我看你是按照书上的步骤一步一步做的，其实有些没必要，建议你看一下傅志方的模态分析第一章。。然后是newmark法，我不可能去一步步给你检查公式的，可以告诉你国内有些著作中公式写错了。。给你推荐（introduction to infiniti element analsys of vibration） dw04116 发表于 2013-7-10 09:32 static/image/common/back.gif
很难一一去看。
我感觉的你搞的太复杂了，就算是想用有限元法也没有必要去计算单元矩阵，直接套用现有的 ...
你可以用你的程序给我讲个例子吗？拜托了 ！
我是用的现有的单元矩阵啊，只是进行了系统的组装而已
页:
[1]
2