论坛 › J01：有限元数值方法 › 外力作用下弹性杆一维应力波传播问题

whutzll 发表于 2008-3-28 10:42:50

外力作用下弹性杆一维应力波传播问题

本程序对照《应力波基础（第二版）》王礼立编，国防工业出版社，北京，2005.8 P339~346
相关公式和说明编写，但是根据结果判定此程序是错误的。
本人才疏学浅，希望得到大家的帮助，先谢谢您的光顾了！
下面有几个疑问，本人认为本程序中有所不足：
1.外力没做处理成节点荷载，只是在{Q}=transpose
2.在固定端没做约束，处理方法如同1
3.结果不正确，平均荷载出入更大
4.有限差分法公式编写是否正确？
5.刚度矩阵和质量矩阵按照书上给的矩阵是否合适？

------------------------------------------程序部分：---------------------------------------------------------------------------
%关键词说明：
%      gNode -------- 节点定义
%      gElement ----- 单元定义
%      gMaterial ---- 材料定义，包括弹性模量，梁的截面积
%      gBC1 --------- 约束条件
%      gDeltaT ------ 时间步长
%      gTimeEnd ----- 计算结束时刻
%      gDisp -------- 位移时程响应
%      gVelo -------- 速度时程响应
%      gAcce -------- 加速度时程响应

function wave_1
% 定义平面杆系的有限元模型数据：

    PlaneFrameModel ;             % 定义有限元模型
    SolveModel ;                  % 求解有限元模型
    SaveResults('wave_1.mat') ;% 保存计算结果
    DisplayResults ;            % 显示计算结果
return ;

function PlaneFrameModel
%定义平面杆系的有限元模型
    global gNode gElement gMaterial gBC1 gDeltaT gTimeEnd gDisp gVelo gAcce
    %单元数目
    n = 3 ;            
    % 节点坐标
    gNode = [ 0      0
             10      0
             20      0
             30      0];

    % 单元定义
    gElement = zeros( n, 3 ) ;
    for i=1:n
      gElement( i, : ) = [ i, i+1, 1 ] ;
    end
   
    % 材料性质
    %         弹性模量    截面积   密度
    gMaterial = ;   %材料 1

    % 第一类约束条件
    %   节点号   自由度号    约束值
    gBC1 = [ 1,      1,      0.0] ;
   
    gDeltaT = 0.001 ;
    gTimeEnd = 3*gDeltaT;    % 计算时间
    timestep = floor(gTimeEnd/gDeltaT) ;

    % 定义位移，速度和加速度
    = size( gNode ) ;
    gDisp = zeros( node_number, timestep ) ;
    gVelo = zeros( node_number, timestep ) ;
    gAcce = zeros( node_number, timestep ) ;
   
    % 初始条件
    gDisp(:,1) = zeros( node_number, 1 ) ;
    gVelo(:,1) = zeros( node_number, 1 ) ;
    gAcce(:,1) = zeros( node_number, 1 ) ;
return

function SolveModel
%求解有限元模型
    global gNode gElement gMaterial gBC1 gK gM gDeltaT gTimeEnd gDisp gVelo gAcce

    % step1. 定义整体刚度矩阵和节点力向量
    = size( gNode ) ;
    gK = sparse( node_number, node_number ) ;
    gM = sparse( node_number, node_number ) ;

    % step2. 计算单元刚度和质量矩阵，并集成到整体刚度和质量矩阵中
    = size( gElement ) ;
    for ie=1:1:element_number
      k = StiffnessMatrix( ie ) ;
      m = MassMatrix( ie );
      AssembleGlobalMatrix( ie, k, m );
    end

    % step3. 计算时程响应
    % step3.1 初始计算
    f0 = sparse( node_number, 1) ;
    f0( node_number , 1 ) = 1.0 ;
    = size( gK ) ;
    K1 = gK + 6/(gDeltaT^2)*gM ;
    timestep = floor(gTimeEnd/gDeltaT) ;
    K1
    % 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)*2 + d ;
    %    K1im(:,ibc) = K1(:,m) ;
    %    K1(:,m) = zeros( node_number, 1 ) ;
    %   K1(m,：） = zeros( 1, node_number ) ;
    %   K1(m,m) = 1.0 ;
    % end
   
    = lu(K1) ;   % 进行三角分解解时间

    % step3.3 对每一个时间步计算
    for i=2:1:timestep
      f1 = f0 + 6/(gDeltaT^2)*gM*( gDisp(:,i-1)+gDeltaT*gVelo(:,i-1)...
                + (gDeltaT^2)/3*gAcce(:,i-1) ) ;
         f1
      % 对f1进行边界条件处理
      % = size( gBC1 ) ;
      %for ibc=1:1:bc1_number
      %   n = gBC1(ibc, 1 ) ;
      %      d = gBC1(ibc, 2 ) ;
      %      m = (n-1)*2 + d ;
      %      f1 = f1 - gBC1(ibc,3) * K1im(:,ibc) ;
      %      f1(m) = gBC1(ibc,3) ;
      %end
      y = KL\f1 ;
      gDisp(:,i) = KU\y ;
      gVelo(:,i) = 3/gDeltaT*(gDisp(:,i)-gDisp(:,i-1))-2*gVelo(:,i-1)...
                     -gDeltaT/2*gAcce(:,i-1);
      gAcce(:,i) = 6/(gDeltaT^2)*(gDisp(:,i)-gDisp(:,i-1))-6/gDeltaT...
                     *gVelo(:,i-1)-2*gAcce(:,i-1);
    end
return

function k = StiffnessMatrix( ie )
%计算单元刚度矩阵
    global gNode gElement gMaterial
    k = zeros( 2, 2 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    A = gMaterial( gElement(ie, 3), 2 ) ;
    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*[1   -1 ;
                -1    1 ] ;
return

function m = MassMatrix( ie )
%计算单元质量矩阵
    global gNode gElement gMaterial
    m = zeros( 2, 2 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    A = gMaterial( gElement(ie, 3), 2 ) ;
    ro = 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) ;
    m = ro*A*L/2*[1,0;
                  0,1] ;
return

function AssembleGlobalMatrix( ie, ke, me )
%把单元刚度和质量矩阵集成到整体刚度矩阵
    global gElement gK gM
    for i=1:1:2
      for j=1:1:2
            m = i ;
            n = j ;
            M = gElement(ie,i) ;
            N = gElement(ie,j) ;
            gK(M,N) = gK(M,N) + ke(m,n) ;
            gM(M,N) = gM(M,N) + me(m,n) ;
      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
   
    fprintf( '节点速度\n' ) ;
    fprintf( '节点号    速度    速度   速度\n' ) ;
    = size( gNode ) ;
    for i=1:node_number
      fprintf('%6d      %16.8e    %16.8e   %16.8e\n',...
                  i,   gVelo(i,1), gVelo(i,2),gVelo(i,3) ) ;
    end
return

[ 本帖最后由 whutzll 于 2008-3-28 10:53 编辑 ]

wxc2000 发表于 2010-10-24 21:16:14

也遇到了相同的问题，lamb的问题，波速似乎太快了

bbssbb 发表于 2010-10-26 09:10:06

是不是单元数太少，结构太刚导致的？

查看完整版本: 外力作用下弹性杆一维应力波传播问题