exam6_2.m

这个压缩包是浙江大学徐荣桥主编的<<结构分析有限元与Matlab程序设计>>一书各章的源代码！

function exam6_2
% 本程序为第六章的第二个算例，采用20结点六面体等参元计算端部受横向力的悬臂梁
%  输入参数： 
%       无

    file_in = 'exam6_2.dat' ;
    
    % 检查文件是否存在
    if exist( file_in ) == 0
        disp( sprintf( '错误：文件 %s 不存在', file_in ) )
        disp( sprintf( '程序终止' ) )
        return ;
    end
    
    % 根据输入文件名称生成输出文件名称
    [path_str,name_str,ext_str] = fileparts( file_in ) ;
    ext_str_out = '.mat' ;
    file_out = fullfile( path_str, [name_str, ext_str_out] ) ;
    
    % 检查输出文件是否存在
    if exist( file_out ) ~= 0 
        answer = input( sprintf( '文件 %s 已经存在，是否覆盖? ( Yes / [No] ):  ', file_out ), 's' ) ;
        if length( answer ) == 0
            answer = 'no' ;
        end
        
        answer = lower( answer ) ;
        if answer ~= 'y' | answer ~= 'yes' 
            disp( sprintf( '请使用另外的文件名，或备份已有的文件' ) ) ;
            disp( sprintf( '程序终止' ) );
            return ; 
        end
    end

    FemModel( file_in ) ;       % 定义有限元模型
    SolveModel ;                % 求解有限元模型
    SaveResults( file_out ) ;   % 保存计算结果
    DisplayResults ;            % 把计算结果显示出来 
return

function FemModel( file_in )
%  定义有限元模型
%  输入参数：
%      file_in ---- 有限模型数据文件
%  返回值：
%      无

%  全局变量及名称
%      gNode ------ 节点坐标
%      gElement --- 单元定义
%      gMaterial -- 材料性质
%      gBC1 ------- 第一类约束条件
%      gNF -------- 集中力
    global gNode gElement gMaterial gBC1 gNF
    
    % 打开文件
    fid = fopen( file_in, 'r' ) ;
    
    % 节点坐标
    node_number = fscanf( fid, '%d', 1 ) ;
    gNode = zeros( node_number, 3 ) ;
    for i = 1:node_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gNode( i, : ) = fscanf( fid, '%f', [1 3] ) ;
    end
    
    % 单元定义
    element_number = fscanf( fid, '%d', 1 ) ;
    gElement = zeros( element_number, 21 ) ;
    for i = 1:element_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gElement( i, : ) = fscanf( fid, '%d', [1 21] ) ;
    end
    
    % 材料性质
    material_number = fscanf( fid, '%d', 1 ) ;
    gMaterial = zeros( material_number, 2 ) ;
    for i = 1:material_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gMaterial( i, : ) = fscanf( fid, '%f', [1 2] ) ;
    end
    
    % 第一类约束条件
    bc1_number = fscanf( fid, '%d', 1 ) ;
    gBC1 = zeros( bc1_number, 3 ) ;
    for i=1:bc1_number
        gBC1(i,1) = fscanf( fid, '%d', 1 ) ;  % 结点号
        gBC1(i,2) = fscanf( fid, '%d', 1 ) ;  % 自由度号
        gBC1(i,3) = fscanf( fid, '%f', 1 ) ;  % 约束值
    end
    
    % 集中力
    nf_number = fscanf( fid, '%d', 1 ) ;
    gNF = zeros( nf_number, 3 ) ;
    for i=1:nf_number
        gNF(i,1) = fscanf( fid, '%d', 1 ) ;  % 结点号
        gNF(i,2) = fscanf( fid, '%d', 1 ) ;  % 自由度号
        gNF(i,3) = fscanf( fid, '%f', 1 ) ;  % 集中力大小
    end
    
    % 关闭文件
    fclose( fid ) ;
return

function SolveModel
%  求解有限元模型
%  输入参数：
%     无
%  返回值：
%     无
%  说明：
%      该函数求解有限元模型，过程如下
%        1. 计算单元刚度矩阵，集成整体刚度矩阵
%        2. 计算单元的等效节点力，集成整体节点力向量
%        3. 处理约束条件，修改整体刚度矩阵和节点力向量
%        4. 求解方程组，得到整体节点位移向量

    global gNode gElement gMaterial gBC1 gNF gK gDelta gElementStress

    % step1. 定义整体刚度矩阵和节点力向量
    [node_number,dummy] = size( gNode ) ;
    gK = sparse( node_number * 3, node_number * 3 ) ;
    f = sparse( node_number * 3, 1 ) ;    

    % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
    [element_number,dummy] = size( gElement ) ;
    for ie=1:1:element_number
        disp( sprintf(  '计算单元刚度矩阵并集成，当前单元: %d', ie  ) ) ;
        k = StiffnessMatrix( ie ) ;
        AssembleStiffnessMatrix( ie, k ) ;
    end
    
    % step3. 把集中力集成到整体节点力向量中
    [nf_number, dummy] = size( gNF ) ;
    for idf = 1:1:nf_number
        node  = gNF(idf, 1) ;
        dim   = gNF(idf, 2) ;
        force = gNF(idf, 3) ;
        f( (node-1)*3+dim ) = f( (node-1)*3+dim ) + force ;
    end
    
    % step4. 处理第一类约束条件，修改刚度矩阵和节点力向量。采用乘大数法
    [bc1_number,dummy] = size( gBC1 ) ;
    for ibc=1:1:bc1_number
        n = gBC1(ibc, 1 ) ;
        d = gBC1(ibc, 2 ) ;
        m = (n-1)*3 + d ;
        f(m,:) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
        gK(m,m) = gK(m,m) * 1e15 ;
    end
    
    % step5. 求解方程组，得到节点位移向量
    gDelta = gK \ f ;
    
    % step6. 计算每个单元的结点应力
    gElementStress = zeros( element_number, 20, 6 ) ;
    delta = zeros( 60, 1 ) ;
    x = [-1, 1, 1, -1, -1, 1, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, -1, 1, 1, -1];
    e = [-1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 0, 0, 0, 0, -1, -1, 1, 1];
    z = [-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 0, 0, 0, 0];
    for ie = 1:element_number
        fprintf( '计算单元应力，当前单元号：%d\n', ie ) ;
        for n=1:20
            B = MatrixB( ie, x(n), e(n), z(n) ) ;
            D = MatrixD( ie ) ;
            delta(1:3:58) = gDelta( gElement(ie,1:20)*3-2) ;
            delta(2:3:59) = gDelta( gElement(ie,1:20)*3-1) ;
            delta(3:3:60) = gDelta( gElement(ie,1:20)*3) ;
            sigma = D*B*delta ;
            gElementStress( ie, n, : ) = sigma ;
        end
    end
return

function k = StiffnessMatrix( ie )   
%  计算平面应变等参数单元的刚度矩阵
%  输入参数：
%     ie -- 单元号
%  返回值：
%     k  -- 单元刚度矩阵
%  说明：
%     用高斯积分法求解平面等参数单元的刚度矩阵

    k = zeros( 60, 60 ) ;
    D = MatrixD( ie ) ;
    %3 x 3 x 3 高斯积分点和权系数
    %x = [-0.774596669241483,              0.0,0.774596669241483] ;   
    %w = [ 0.555555555555556,0.888888888888889,0.555555555555556] ;  
    %for i=1:1:length(x)
    %    for j=1:1:length(x)
    %        for m=1:1:length(x)
    %            B = MatrixB( ie, x(i), x(j), x(m) ) ;
    %            J = Jacobi( ie, x(i), x(j), x(m) ) ;
    %            k = k + w(i)*w(j)*w(m)*transpose(B)*D*B*det(J) ;   
    %        end
    %    end
    %end
    
    %14点 irons 积分, 详细参见一下文献 
    % Irons B. M., 1971 Quadrature rules for brick based finite elements, 
    %   International Journal for Numerical Methods in Engineering, 3(2):293-294.
    b0 = 320/361 ;
    c0 = 121/361 ;
    b  = sqrt( 19/30 ) ;
    c  = sqrt( 19/33 ) ;
    x1 = [-b,  0,  0
           b,  0,  0
           0, -b,  0
           0,  b,  0
           0,  0, -b
           0,  0,  b ] ;
    x2 = [-c, -c, -c
           c,  c,  c
           c, -c, -c
          -c,  c, -c
          -c, -c,  c
           c,  c, -c
          -c,  c,  c
           c, -c,  c ] ;
    for i=1:6
        B = MatrixB( ie, x1(i,1), x1(i,2), x1(i,3) ) ;
        J = Jacobi( ie, x1(i,1), x1(i,2), x1(i,3) ) ;
        k = k + b0*transpose(B)*D*B*det(J) ;
    end
    for i=1:8
        B = MatrixB( ie, x2(i,1), x2(i,2), x2(i,3) ) ;
        J = Jacobi( ie, x2(i,1), x2(i,2), x2(i,3) ) ;
        k = k + c0*transpose(B)*D*B*det(J) ;
    end
return

function D = MatrixD( ie )
%  计算单元的弹性矩阵D
%  输入参数：
%     ie --------- 单元号
%  返回值：
%     D  --------- 弹性矩阵D

    global gElement gMaterial 
    E  = gMaterial( gElement(ie, 21), 1 ) ;   % 弹性模量
    mu = gMaterial( gElement(ie, 21), 2 ) ;   % 泊松比 
    A1 = mu/(1-mu) ;                        
    A2 = (1-2*mu)/2/(1-mu) ;                
    A3 = E*(1-mu)/(1+mu)/(1-2*mu) ;         
    D = A3*[  1  A1  A1   0   0   0
             A1   1  A1   0   0   0
             A1  A1   1   0   0   0
              0   0   0  A2   0   0
              0   0   0   0  A2   0
              0   0   0   0   0  A2] ;
return

function B = MatrixB( ie, xi, eta, zeta )
%  计算单元的应变矩阵B
%  输入参数：
%     ie --------- 单元号
%     xi,eta ----- 局部坐标  
%  返回值：
%     B  --------- 在局部坐标处的应变矩阵B

    [N_x, N_y, N_z] = N_xyz( ie, xi, eta, zeta ); 
    B = zeros( 6, 60 ) ;
    for i=1:1:20
        B(:,(3*i-2):3*i) = [ N_x(i)       0       0
                                  0  N_y(i)       0
                                  0       0  N_z(i)
                                  0  N_z(i)  N_y(i)
                             N_z(i)      0   N_x(i)
                             N_y(i)  N_x(i)       0 ] ;
    end