本帖最后由 ME! 于 2013-7-8 21:48 编辑
下面是我的程序，能帮我看看吗，共振处的频率不等于omega1，我用振型叠加法也是得不到共振频率，不知道是哪里错了- %-------------------------------------------------------------------------%
- %--------------------------------------------------------------------------
- %newmark
- % (0) 输入控制数据
- %--------------------------------------------------------------------------
- clear; clc;
- jk=3;
- c= [4,13]; %约束的节点 两点支撑，4为前轴承为节点，13为后轴承位节点
- jth=14; % 第jth自由度方向的响应
- [n1 n2]=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= [16.4 12.6 10 10 6.5 3 6.5 47.2 5 5 5 10.5 10 10 6.3 23 ]; %长度分段
- d1= [15 16 17 17 20.5 29 21.7 21 20 20.5 20.5 20.5 17 17 16 12.5]; %主轴分段外径
- d2= [12 9 9 8 8 8 8 8 8 8 8 8 8 8 8 8 ]; %主轴分段内径
- dm1=[0 24 35 35 0 0 0 30 0 0 0 0 35 35 28 25 ]; %主轴分段附加零件外径
- dm2=[0 16 17 17 0 0 0 21 0 17 17 0 17 17 16 12.5]; %主轴分段附加零件内径
-
- 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; %中间轴承的径向刚度
- [KK1,KK2,KK3]=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;%将后轴承径向结合部刚度值并入整体刚度矩阵
- [bc1_number,dummy] = 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
- [gEigVector, gEigValue] = 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=[100,1,2]; % 施加激振力的幅值
- 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 ) ) ;
- [N,N] = 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进行边界条件处理
- [bc1_number,dummy] = 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
- [KL,KU] = 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进行边界条件处理
- [bc1_number,dummy] = 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,:)
- %--------------------------------------------------------------------------
- % 绘制时程曲线
- %[node_number,dummy] = 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',[1 1 1]);
- %set(0,'defaultaxesfontsize',8)
- xlabel( '频率(Hz)') ;
- ylabel( '幅度' ) ;
-
|
发表于 2012-6-1 11:39
提升卡
变色卡
显身卡
楼主
发表于 2012-6-1 20:23
