ME! 发表于 2012-12-14 21:52

给大家贴一个均匀梁单元的matlab求振型和响应曲线的程序

function exam8_1
% 本程序为第八章的第一个算例,采用平面梁单元计算两铰抛物线拱的自由振动特性
%      输入参数: 无
%      输出结果: 前3阶振动频率及其相应的振型
% 定义全局变量
%      gNode ------ 节点坐标
%      gElement --- 单元定义
%      gMaterial -- 材料性质
%      gBC1 ------- 第一类约束条件
%      gK --------- 整体刚度矩阵
%      gDelta ----- 整体节点坐标
    PlaneFrameModel ;      % 定义有限元模型
    SolveModel ;         % 求解有限元模型
    DisplayResults ;       % 显示计算结果
return ;
function PlaneFrameModel
%定义平面杆系的有限元模型
%输入参数:
%      无
%返回值:
%      无
%说明:
%      该函数定义平面杆系的有限元模型数据:
%      gNode ------- 节点定义
%      gElement ---- 单元定义
%      gMaterial --- 材料定义,包括弹性模量,梁的截面积和梁的抗弯惯性矩
%      gBC --------- 约束条件
    global gNode gElement gMaterial gBC1
    % 给定抛物线拱的几何特征
    L = 60 ;               %计算跨径   
    f = 7.5 ;            %计算矢高
   
    n = 100 ;            %单元数目
    x = -L/2:L/n:L/2 ;   %结点的x坐标
    a = f/L^2*4 ;
    y = - a * x.^2 ;       %结点的y坐标
    % 节点坐标
    gNode = ;
   
    % 单元定义
    gElement = zeros( n, 3 ) ;
    for i=1:n
      gElement( i, : ) = [ i, i+1, 1 ] ;
    end
   
    % 材料性质
    %         弹性模量   抗弯惯性矩   截面积   密度
    gMaterial = ;   %材料 1
    % 第一类约束条件
    %   节点号   自由度号    约束值
    gBC1 = [ 1,      1,      0.0
             1,      2,      0.0
             n+1,      1,      0.0
             n+1,      2,      0.0] ;
return
function SolveModel
%求解有限元模型
%输入参数:
%   无
%返回值:
%   无
%说明:
%      该函数求解有限元模型,过程如下
%      1. 计算单元的刚度和质量矩阵,集成整体刚度和质量矩阵
%      2. 处理约束条件,修改整体刚度矩阵
%      3. 求解特征值问题
    global gNode gElement gMaterial gBC1 gK gM gEigValue gEigVector
    % step1. 定义整体刚度矩阵和节点力向量
    %sparse创建稀疏矩阵S = sparse(i,j,s,m,n,nzmax)
    %由向量i,j,s生成一个m*n的含有nzmax个非零元素的稀疏矩阵S,并且有 S(i(k),j(k)) = s(k)。向量 i,j 和 s 有相同的长度。
    %S = sparse(i,j,s,m,n)这种情况下nzmax = length(s)。
    gM = sparse( node_number * 3, node_number * 3 ) ;
    = size( gNode ) ;
    gK = sparse( node_number * 3, node_number * 3 ) ;   

    % step2. 计算单元刚度和质量矩阵,并集成到整体刚度和质量矩阵中
    = size( gElement ) ;
    for ie=1:1:element_number
      k = StiffnessMatrix( ie ) ;
      m = MassMatrix( ie ) ;
      AssembleGlobalMatrix( ie, k, m ) ;
    end
    % step3. 处理第一类约束条件,修改刚度矩阵和质量矩阵。(采用划行划列法)
    = size( gBC1 ) ;
    w2max = max( diag(gK)./diag(gM) ) ;
    for ibc=1:1:bc1_number
      n = gBC1(ibc, 1 ) ;
      d = gBC1(ibc, 2 ) ;
      m = (n-1)*3 + d ;
      gK(:,m) = zeros( node_number*3, 1 ) ;
      gK(m,:) = zeros( 1, node_number*3 ) ;
      gK(m,m) = 1;
      gM(:,m) = zeros( node_number*3, 1 ) ;
      gM(m,:) = zeros( 1, node_number*3 ) ;
      gM(m,m) = gK(m,m)/w2max/1e10 ;
    end
   
    % step4. 求解特征值问题
    % step4.1为了使刚度矩阵和质量矩阵对称(在计算时可能引入舍入误差)
    for i=1:node_number*3
      for j=i:node_number*3
            gK(j,i) = gK(i,j) ;
            gM(j,i) = gM(i,j) ;
      end
    end
   
    % step4.2 计算前6阶特征值和特征向量
    = eigs(gK, gM, 3, 'SM' ) ;
   
    % step4.3 修改特征向量中受约束的自由度
    for ibc=1:1:bc1_number
      n = gBC1(ibc, 1 ) ;
      d = gBC1(ibc, 2 ) ;
      m = (n-1)*3 + d ;
      gEigVector(m,:) = gBC1(ibc,3) ;
    end
return
function k = StiffnessMatrix( ie )
%计算单元刚度矩阵
%输入参数:
%   ie -------单元号
%返回值:
%   k----整体坐标系下的刚度矩阵
    global gNode gElement gMaterial
    k = zeros( 6, 6 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    I = gMaterial( gElement(ie, 3), 2 ) ;
    A = gMaterial( gElement(ie, 3), 3 ) ;
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
    k = [E*A/L         0          0 -E*A/L         0          0
               012*E*I/L^36*E*I/L^2      0 -12*E*I/L^36*E*I/L^2
               0   6*E*I/L^2    4*E*I/L      0-6*E*I/L^2    2*E*I/L
          -E*A/L         0          0E*A/L         0          0
               0 -12*E*I/L^3 -6*E*I/L^2      012*E*I/L^3 -6*E*I/L^2
               0   6*E*I/L^2    2*E*I/L      0-6*E*I/L^2    4*E*I/L] ;
    T = TransformMatrix( ie ) ;
    k = T*k*transpose(T) ;
return
function m = MassMatrix( ie )
%计算单元质量矩阵
%输入参数:
%   ie -------单元号
%返回值:
%   m----整体坐标系下的质量矩阵
    global gNode gElement gMaterial
    m = zeros( 6, 6 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    A = gMaterial( gElement(ie, 3), 3 ) ;
    ro = gMaterial( gElement(ie, 3 ), 4 ) ;
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
    m = ro*A*L/420*[140      0      0   70      0      0
                      0    156   22*L    0   54   -13*L
                      0   22*L4*L^2    0   13*L-3*L^2
                     70      0      0140      0       0
                      0   54   13*L    0    156   -22*L
                      0-13*L   -3*L    0-22*L4*L^2 ] ;
    T = TransformMatrix( ie ) ;
    m = T*m*transpose(T) ;
return
function AssembleGlobalMatrix( ie, ke, me )
%把单元刚度和质量矩阵集成到整体刚度矩阵
%输入参数:
%      ie--- 单元号
%      ke--- 单元刚度矩阵
%      me--- 单元质量矩阵
%返回值:
%      无
    global gElement gK gM
    for i=1:1:2
      for j=1:1:2
            for p=1:1:3
                for q =1:1:3
                  m = (i-1)*3+p ;
                  n = (j-1)*3+q ;
                  M = (gElement(ie,i)-1)*3+p ;
                  N = (gElement(ie,j)-1)*3+q ;
                  gK(M,N) = gK(M,N) + ke(m,n) ;
                  gM(M,N) = gM(M,N) + me(m,n) ;
                end
            end
      end
    end
return
function T = TransformMatrix( ie )
%计算单元的坐标转换矩阵( 局部坐标 -> 整体坐标 )
%输入参数
%      ie----- 节点号
%返回值
%      T ------- 从局部坐标到整体坐标的坐标转换矩阵
    global gElement gNode
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = sqrt( (xj-xi)^2 + (yj-yi)^2 ) ;
    c = (xj-xi)/L ;
    s = (yj-yi)/L ;
    T=[ c-s   0   0   0   0
      s   c   0   0   0   0
      0   0   1   0   0   0
      0   0   0   c-s   0
      0   0   0   s   c   0
      0   0   0   0   0   1] ;
return
function DisplayResults
%显示计算结果
%输入参数:
%   无
%返回值:
%   无
    global gNode gElement gMaterial gBC1 gEigValue gEigVector
    fre_number = length(diag(gEigValue)) ;
   
    % 打印特征向量(振型)
    fprintf( '\n\n 表一   特征向量(振型)\n' ) ;
    for i=1:fre_number
      fprintf( '----------------') ;
    end
    fprintf( '\n' ) ;
    for i=1:fre_number
      fprintf( '%6d      ', i ) ;
    end
    fprintf( '\n' ) ;
    for i=1:fre_number
      fprintf( '----------------') ;
    end
    fprintf( '\n' ) ;
    =size(gEigVector) ;
    for i=1:dof
      for j=fre_number:-1:1
            fprintf( '%15.7e ', gEigVector(i,j) ) ;
      end
      fprintf( '\n' ) ;
    end
    for i=1:fre_number
      fprintf( '----------------') ;
    end
    fprintf( '\n' ) ;
    % 打印特征值
    fprintf( '\n\n\n\n 表二   特征值(频率)列表\n' ) ;
    fprintf( '----------------------------------------------------------\n') ;
    fprintf( '   阶数      特征值          频率(Hz)      圆频率(Hz)\n' ) ;
    fprintf( '----------------------------------------------------------\n') ;
    for i=fre_number:-1:1
      fprintf( '%6d   %15.7e   %15.7e   %15.7e\n', fre_number-i+1, ...
            gEigValue(i,i), sqrt(gEigValue(i,i))/2/pi, sqrt(gEigValue(i,i)) ) ;
    end
    fprintf( '----------------------------------------------------------\n') ;
   
    % 绘制振型图
    for j=fre_number:-1:1
      figure ;
      x = gNode(:,1) ;
      y = gNode(:,2) ;
      dx = gEigVector(1:3:length(x)*3, j ) ;
      dy = gEigVector(2:3:length(x)*3, j ) ;
      factor = max( )* 0.05;
      plot(x,y,'-', x+factor*dx, y+factor*dy, ':') ;
      title( sprintf( '第%d阶频率: %.3f Hz', fre_number-j+1, sqrt(gEigValue(j,j))/2/pi ) ) ;
      axis equal;
      axis off ;
    end
return






function exam8_2
% 本程序为第八章的第二个算例,采用平面梁单元计算两铰抛物线拱的在初始条件下
%自由振动,并对时程曲线结果进行FFT变换,求得的频率可与exam8_1.m的结果进行
%比较,以验证本程序的可靠性
%      输入参数: 无
%      输出结果: 位移,速度和加速度的时程曲线
% 定义全局变量
%      gNode ------ 节点坐标
%      gElement --- 单元定义
%      gMaterial -- 材料性质
%      gBC1 ------- 第一类约束条件
%      gK --------- 整体刚度矩阵
%      gDelta ----- 整体节点坐标
    PlaneFrameModel ;             % 定义有限元模型
    SolveModel ;                  % 求解有限元模型
    SaveResults('exam8_2.mat') ;% 保存计算结果
    DisplayResults ;            % 显示计算结果
return ;
function PlaneFrameModel
%定义平面杆系的有限元模型
%输入参数:
%      无
%返回值:
%      无
%说明:
%      该函数定义平面杆系的有限元模型数据:
%      gNode -------- 节点定义
%      gElement ----- 单元定义
%      gMaterial ---- 材料定义,包括弹性模量,梁的截面积和梁的抗弯惯性矩
%      gBC1 --------- 约束条件
%      gDeltaT ------ 时间步长
%      gTimeEnd ----- 计算结束时刻
%      gDisp -------- 位移时程响应
%      gVelo -------- 速度时程响应
%      gAcce -------- 加速度时程响应
    global gNode gElement gMaterial gBC1 gDeltaT gTimeEnd gDisp gVelo gAcce
    % 给定抛物线拱的几何特征
    L = 60 ;               %计算跨径(m)   
    f = 7.5 ;            %计算矢高(m)
   
    n = 100 ;            %单元数目
    x = -L/2:L/n:L/2 ;   %结点的x坐标
    a = f/L^2*4 ;
    y = - a * x.^2 ;       %结点的y坐标
    % 节点坐标
    gNode = ;
   
    % 单元定义
    gElement = zeros( n, 3 ) ;
    for i=1:n
      gElement( i, : ) = [ i, i+1, 1 ] ;
    end
   
    % 材料性质
    %         弹性模量   抗弯惯性矩   截面积   密度
    gMaterial = ;   %材料 1
    % 第一类约束条件
    %   节点号   自由度号    约束值
    gBC1 = [ 1,      1,      0.0
             1,      2,      0.0
             n+1,      1,      0.0
             n+1,      2,      0.0] ;
   
    gDeltaT = 0.01 ;
    gTimeEnd = 4096*gDeltaT;    % 计算时间为载荷通过所需时间的两倍
    timestep = floor(gTimeEnd/gDeltaT) ;
    % 定义位移,速度和加速度
    gDisp = zeros( (n+1)*3, timestep ) ;
    gVelo = zeros( (n+1)*3, timestep ) ;
    gAcce = zeros( (n+1)*3, timestep ) ;
   
    % 初始条件
    gDisp(:,1) = zeros( (n+1)*3, 1 ) ;
    gVelo(:,1) = ones( (n+1)*3, 1 ) ;
return
function SolveModel
%求解有限元模型
%输入参数:
%   无
%返回值:
%   无
%说明:
%      该函数求解有限元模型,过程如下
%      1. 计算单元的刚度和质量矩阵,集成整体刚度和质量矩阵
%      2. 用Newmark法计算时程响应
    global gNode gElement gMaterial gBC1 gK gM gDeltaT gTimeEnd gDisp gVelo gAcce
    % step1. 定义整体刚度矩阵和节点力向量
    = size( gNode ) ;
    gK = sparse( node_number * 3, node_number * 3 ) ;
    gM = sparse( node_number * 3, node_number * 3 ) ;
    % step2. 计算单元刚度和质量矩阵,并集成到整体刚度和质量矩阵中
    = size( gElement ) ;
    for ie=1:1:element_number
      k = StiffnessMatrix( ie ) ;
      m = MassMatrix( ie ) ;
      AssembleGlobalMatrix( ie, k, m ) ;
    end
    % step3. 计算时程响应(Newmark法)
    % step3.1 初始计算
    gama = 0.5 ;
    beta = 0.25 ;
    C = zeros( size( gK ) ) ;
    = size( gK ) ;
    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 = gK + alpha0*gM + alpha1*C;
    timestep = floor(gTimeEnd/gDeltaT) ;
   
    % step3.2 对K1进行边界条件处理
    = size( gBC1 ) ;
    K1im = zeros(N,bc1_number) ;
    for ibc=1:1:bc1_number
      n = gBC1(ibc, 1 ) ;
      d = gBC1(ibc, 2 ) ;
      m = (n-1)*3 + d ;
      K1im(:,ibc) = K1(:,m) ;
      K1(:,m) = zeros( node_number*3, 1 ) ;
      K1(m,:) = zeros( 1, node_number*3 ) ;
      K1(m,m) = 1.0 ;
    end
    = lu(K1) ;   % 进行三角分解,节省后面的求解时间
   
    % step3.3 计算初始加速度
    gAcce(:,1) = gM\(-gK*gDisp(:,1)-C*gVelo(:,1)) ;
   
    % step3.4 对每一个时间步计算
    for i=2:1:timestep
      fprintf( '当前时间步:%d\n', i ) ;
      f1 = gM*(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( gBC1 ) ;
      for ibc=1:1:bc1_number
            n = gBC1(ibc, 1 ) ;
            d = gBC1(ibc, 2 ) ;
            m = (n-1)*3 + d ;
            f1 = f1 - gBC1(ibc,3) * K1im(:,ibc) ;
            f1(m) = gBC1(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
return
function k = StiffnessMatrix( ie )
%计算单元刚度矩阵
%输入参数:
%   ie -------单元号
%返回值:
%   k----整体坐标系下的刚度矩阵
    global gNode gElement gMaterial
    k = zeros( 6, 6 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    I = gMaterial( gElement(ie, 3), 2 ) ;
    A = gMaterial( gElement(ie, 3), 3 ) ;
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
    k = [E*A/L         0          0 -E*A/L         0          0
               012*E*I/L^36*E*I/L^2      0 -12*E*I/L^36*E*I/L^2
               0   6*E*I/L^2    4*E*I/L      0-6*E*I/L^2    2*E*I/L
          -E*A/L         0          0E*A/L         0          0
               0 -12*E*I/L^3 -6*E*I/L^2      012*E*I/L^3 -6*E*I/L^2
               0   6*E*I/L^2    2*E*I/L      0-6*E*I/L^2    4*E*I/L] ;
    T = TransformMatrix( ie ) ;
    k = T*k*transpose(T) ;
return
function m = MassMatrix( ie )
%计算单元质量矩阵
%输入参数:
%   ie -------单元号
%返回值:
%   m----整体坐标系下的质量矩阵
    global gNode gElement gMaterial
    m = zeros( 6, 6 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    A = gMaterial( gElement(ie, 3), 3 ) ;
    ro = gMaterial( gElement(ie, 3 ), 4 ) ;
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
    m = ro*A*L/420*[140      0      0   70      0      0
                      0    156   22*L    0   54   -13*L
                      0   22*L4*L^2    0   13*L-3*L^2
                     70      0      0140      0       0
                      0   54   13*L    0    156   -22*L
                      0-13*L   -3*L    0-22*L4*L^2 ] ;
    T = TransformMatrix( ie ) ;
    m = T*m*transpose(T) ;
return
function AssembleGlobalMatrix( ie, ke, me )
%把单元刚度和质量矩阵集成到整体刚度矩阵
%输入参数:
%      ie--- 单元号
%      ke--- 单元刚度矩阵
%      me--- 单元质量矩阵
%返回值:
%      无
    global gElement gK gM
    for i=1:1:2
      for j=1:1:2
            for p=1:1:3
                for q =1:1:3
                  m = (i-1)*3+p ;
                  n = (j-1)*3+q ;
                  M = (gElement(ie,i)-1)*3+p ;
                  N = (gElement(ie,j)-1)*3+q ;
                  gK(M,N) = gK(M,N) + ke(m,n) ;
                  gM(M,N) = gM(M,N) + me(m,n) ;
                end
            end
      end
    end
return
function T = TransformMatrix( ie )
%计算单元的坐标转换矩阵( 局部坐标 -> 整体坐标 )
%输入参数
%      ie----- 节点号
%返回值
%      T ------- 从局部坐标到整体坐标的坐标转换矩阵
    global gElement gNode
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = sqrt( (xj-xi)^2 + (yj-yi)^2 ) ;
    c = (xj-xi)/L ;
    s = (yj-yi)/L ;
    T=[ c-s   0   0   0   0
      s   c   0   0   0   0
      0   0   1   0   0   0
      0   0   0   c-s   0
      0   0   0   s   c   0
      0   0   0   0   0   1] ;
return
function ft = NodeForce( t )
% 计算在时刻 t 的结点载荷列阵
% 输入参数
%    t ------ 时刻
% 返回值
%   ft ------ t时刻的结点载荷列阵
    global gNode gElement gLoad gLoadVelo
    Load = gLoad*sin(2*pi*50*t) ;
    %Load = gLoad ;
    node_number = size( gNode, 1 ) ;
    ft = zeros( node_number*3, 1 ) ;
    xt = gNode(1,1) + gLoadVelo * t ;
    element_number = size( gElement, 1) ;
    for ie=1:element_number
      node = gElement( ie, 1:2 ) ;
      x = gNode( node, 1 ) ;
      y = gNode( node, 2 ) ;
      L = sqrt( (x(2)-x(1))^2 + (y(2)-y(1))^2 ) ;
      cosq = (x(2)-x(1))/L ;
      sinq = (y(2)-y(1))/L ;
      N = -Load*sinq ;
      P = -Load*cosq ;
      if x(1) <= xt & x(2) >=xt
            N1 = (x(2)-xt)/(x(2)-x(1)) * N ;
            N2 = (xt-x(1))/(x(2)-x(1)) * N ;
            dx = (xt-x(1))/cosq ;
            P1 = P*(1-3*(dx/L)^2+2*(dx/L)^3) ;
            P2 = P*( +3*(dx/L)^2-2*(dx/L)^3) ;
            M1 = P*dx*(1-2*dx/L+(dx/L)^2) ;
            M2 = P*dx*( -dx/L + (dx/L)^2) ;
            T = TransformMatrix( ie ) ;
            fe = T * ;
            ft( node(1)*3-2:node(1)*3 ) = fe(1:3) ;
            ft( node(2)*3-2:node(2)*3 ) = fe(4:6) ;
            break ;
      end
    end
return
function SaveResults( file_out )
%保存计算结果
%输入参数:
%   无
%返回值:
%   无
    global gNode gElement gMaterial gBC1 gDeltaT gTimeEnd gLoad gLoadVelo gDisp gVelo gAcce
    save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1',...
          'gDeltaT', 'gTimeEnd', 'gLoad', 'gLoadVelo', 'gDisp', 'gVelo', 'gAcce' ) ;
return
function DisplayResults
%显示计算结果
%输入参数:
%   无
%返回值:
%   无
    global gNode gElement gMaterial gBC1 gDisp gVelo gAcce gDeltaT gTimeEnd
    % 绘制时程曲线
    = size(gNode) ;
    t = 0:gDeltaT:gTimeEnd-gDeltaT ;
    d = gDisp((floor(node_number/4)*3)+2,:) ;
    v = gVelo((floor(node_number/4)*3)+2,:) ;
    a = gAcce((floor(node_number/4)*3)+2,:) ;
    tt = 0:gDeltaT/10:gTimeEnd-gDeltaT ;
    dd = spline(t,d,tt) ;%spline表示三次样条插值
    vv = spline(t,v,tt) ;
    aa = spline(t,a,tt) ;
    figure ;
    plot( tt,dd, tt, vv, tt, aa ) ;
    title( 'L/4处挠度时程曲线' ) ;
    xlabel( '时间(s)') ;
    ylabel( '幅度(m)' ) ;
   
    % 对时程曲线进行FFT变换,获取频谱特性
    fd = fft( d ) ;
    figure ;
    df = 1/gTimeEnd ;%FFT变换的频率分辨率
    f = (0:length(d)-1)*df ;
    plot(f,abs(fd)) ;
    set(gca,'xlim',) ;
    title( 'L/4处挠度的频谱图' ) ;
    xlabel( '频率(Hz)') ;
    ylabel( '幅度' ) ;
return
顺便问问大家如果变截面梁,应该怎么做,期待高手的回答

mkhatebm 发表于 2013-3-7 18:11

兄弟,无论等截面还是变截面,都可以通过上述方法求的,只是把节点参数改一下就可以了

ME! 发表于 2013-3-7 18:40

mkhatebm 发表于 2013-3-7 18:11 static/image/common/back.gif
兄弟,无论等截面还是变截面,都可以通过上述方法求的,只是把节点参数改一下就可以了

呵呵,这个问题早就改好了

自由自在 发表于 2013-3-21 10:35

请教一下,如果不是简单的直梁结构,应该怎么得到它的固有频率和振型啊,比图说"L”型、"U"型的悬臂梁结构

ME! 发表于 2013-3-21 14:37

这个没有碰到过,但是我觉得是一样的吧!

yshw2008 发表于 2013-3-31 09:44

想问下你指的是那一本书啊,学习下

ME! 发表于 2013-4-1 09:47

yshw2008 发表于 2013-3-31 09:44 static/image/common/back.gif
想问下你指的是那一本书啊,学习下

?你是说例子吗?

mwglikai 发表于 2013-4-2 09:39

ME! 发表于 2013-4-1 09:47 static/image/common/back.gif
?你是说例子吗?

如果是轴承约束,有四个刚度系数,怎么处理?谢谢!

ME! 发表于 2013-4-2 16:03

mwglikai 发表于 2013-4-2 09:39 static/image/common/back.gif
如果是轴承约束,有四个刚度系数,怎么处理?谢谢!

你是说一个轴承的四个刚度系数,包括耦合刚度吗?因为我用单元矩阵的时候,节点位移我设为0了,所以加上轴承刚度也不会有影响了
还有你是毕业了再湘潭工作吗?我也是湖大的哦!

mgh_nx 发表于 2013-5-5 00:58

自由自在 发表于 2013-3-21 10:35 static/image/common/back.gif
请教一下,如果不是简单的直梁结构,应该怎么得到它的固有频率和振型啊,比图说"L”型、"U"型的悬臂梁结构

横截面刚度改变一下就可以了。

mxlzhenzhu 发表于 2013-5-5 12:33

感谢分享;

变截面梁,呵呵,这就是用这种单元的麻烦的地方,为什么不直接用实体单元?
我走的是捷径,用matlab直接读取nastran生成的矩阵文件,省去了自己组装的过程。

ME! 发表于 2013-5-6 10:13

mxlzhenzhu 发表于 2013-5-5 12:33 static/image/common/back.gif
感谢分享;

变截面梁,呵呵,这就是用这种单元的麻烦的地方,为什么不直接用实体单元?

其实组装也没有那么麻烦!你是说ansys的实体单元吗?自由度不一样,我计算出来结果相差很大,所以不敢用ansys实体单元,但是ansys的beam3单元自由度和平面梁的自由度是一样的,但是我用beam3梁加上轴承弹簧刚度之后对结果没有什么影响,难道这就是它的弊端?

紫域云庭 发表于 2013-10-27 19:54

你好,请问 这两个例子 是什么参考书上的例子,能否发给我啊~
我邮箱ybliang2010@gmail.com谢谢了

ME! 发表于 2013-10-28 19:22

紫域云庭 发表于 2013-10-27 19:54 static/image/common/back.gif
你好,请问 这两个例子 是什么参考书上的例子,能否发给我啊~
我邮箱    谢谢了

徐斌--Matlab有限元结构动力学分析与工程应用-

402144999 发表于 2014-7-26 10:53

页: [1] 2
查看完整版本: 给大家贴一个均匀梁单元的matlab求振型和响应曲线的程序