exam6_1.m

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

function exam6_1
% 本程序为第六章的第一个算例，采用4～8结点四边形等参元计算受内压的旋转厚壁圆筒
%  输入参数： 
%       无

    % 输出文件名称
    file_out = 'exam6_1.mat' ;
    
    % 检查输出文件是否存在
    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 ;                  % 定义有限元模型
    SolveModel ;                % 求解有限元模型
    SaveResults( file_out ) ;   % 保存计算结果
    DisplayResults ;            % 把计算结果显示出来 
return

function FemModel
%  定义有限元模型
%  输入参数：
%      无
%  返回值：
%      无

%  全局变量及名称
%      gNode ------ 节点坐标
%      gElement --- 单元定义
%      gMaterial -- 材料性质
%      gBC1 ------- 第一类约束条件
%      gBC3 ------- 斜约束
%      gDF -------- 分布力
    global gNode gElement gMaterial gBC1 gBC3 gDF
    
    r1 = 0.05 ;              %  内半径
    r2 = 0.1 ;               %  外半径
    theta = 10.0*pi/180.0 ;  %  圆心角（弧度）
    
    % 节点坐标
    node_number = 28 ;
    gNode = zeros( node_number, 2 ) ;
    j = 0 ;
    for i = 1:5:26     % 结点1,6,11,16,21,26的坐标
        r = r1 + (r2-r1)/5 * j ;
        gNode(i,1) = r*cos(theta) ;
        gNode(i,2) = r*sin(theta) ;
        j = j + 1 ;
    end
    j = 0 ;
    for i = 4:5:24     % 结点4,9,14,19,24的坐标
        r = r1+(r2-r1)/10+(r2-r1)/5*j ;
        gNode(i,1) = r*cos(theta) ;
        gNode(i,2) = r*sin(theta) ;
        j = j+1 ;
    end
    j = 0 ;
    for i=2:5:27       % 结点2,7,12,17,22,27的坐标
        r = r1+(r2-r1)/5*j ;
        gNode(i,1) = r*cos(theta/2) ;
        gNode(i,2) = r*sin(theta/2) ;
        j = j+1 ;
    end
    j = 0 ;
    for i=3:5:28      % 结点3,8,13,18,23,28的坐标
        r = r1+(r2-r1)/5*j ;
        gNode(i,1) = r ;
        gNode(i,2) = 0 ;
        j = j+1 ;
    end
    j = 0 ;
    for i=5:5:25      % 结点5,10,15,20,25的坐标
        r = r1+(r2-r1)/10+(r2-r1)/5*j ;
        gNode(i,1) = r ;
        gNode(i,2) = 0 ;
        j = j+1 ;
    end
    
    % 单元定义
    element_number = 5 ;
    gElement = zeros( element_number, 9 ) ;
    gElement(1,:)=[3  8  6  1   5  7  4  2  1] ;
    for i=2:5
        gElement(i,1:8) = gElement(i-1,1:8)+5 ;
        gElement(i,9) = 1 ;
    end
    
    % 材料性质
    gMaterial = [2.0e11 0.3 7800] ;
    
    % 第一类约束条件
    bc1_number = 11;
    gBC1 = zeros( bc1_number, 3 ) ;
    j = 0 ;
    for i=3:5:28      %  结点3,8,13,18,23,28的边界条件(y方向约束)
        j = j + 1 ;
        gBC1(j,1) = i ;  % 结点号
        gBC1(j,2) = 2 ;  % 自由度号
        gBC1(j,3) = 0 ;  % 约束值
    end
    j = 6 ;
    for i=5:5:25     %  结点5,10,15,20,25的边界条件(y方向约束)
        j = j + 1;
        gBC1(j,1) = i ;  % 结点号
        gBC1(j,2) = 2 ;  % 自由度号
        gBC1(j,3) = 0 ;  % 约束值
    end
    
    % 斜约束条件
    bc3_number = 17 ;
    gBC3 = zeros( bc3_number, 4 ) ;
    j = 0 ;
    for i=1:5:26      %  结点1,6,11,16,21,26的边界条件(y'方向约束)
        j = j + 1 ;
        gBC3(j,1) = i ;     % 结点号
        gBC3(j,2) = 2 ;     % 自由度号
        gBC3(j,3) = 0 ;     % 约束值
        gBC3(j,4) = theta ; % 倾斜角度
    end
    j = 6 ;
    for i=4:5:24     %  结点4,9,14,19,24的边界条件(y'方向约束)
        j = j + 1;
        gBC3(j,1) = i ;     % 结点号
        gBC3(j,2) = 2 ;     % 自由度号
        gBC3(j,3) = 0 ;     % 约束值
        gBC3(j,4) = theta ; % 倾斜角度
    end
    j = 11 ;
    for i=2:5:27 
        j = j+1 ;
        gBC3(j,1) = i ;       % 结点号
        gBC3(j,2) = 2 ;       % 自由度号
        gBC3(j,3) = 0 ;       % 约束值
        gBC3(j,4) = theta/2 ; % 倾斜角度
    end
    
    % 分布压力
    df_number = 1 ;
    gDF = [3 1 2 120e6 120e6 120e6] ;
return

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

    global gNode gElement gMaterial gBC1 gBC3 gDF gK gDelta1 gDelta gElementStress

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

    % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
    [element_number,dummy] = size( gElement ) ;
    for ie=1:1:element_number
        disp( sprintf(  '计算单元刚度矩阵并集成，当前单元: %d', ie  ) ) ;
        k = StiffnessMatrix( ie ) ;
        AssembleStiffnessMatrix( ie, k ) ;
    end
    
    % step3. 计算单元离心力的等效结点力，并集成到整体结点力向量中
    for ie=1:1:element_number
        evf = EquivalentVolumeForce( ie ) ;
        node = gElement( ie, 1:8 ) ;
        f( (node-1)*2+1,2 ) = f( (node-1)*2+1,2 ) + evf(1:2:15) ;
        f( (node-1)*2+2,2 ) = f( (node-1)*2+2,2 ) + evf(2:2:16) ;
    end

    % step4. 计算分布压力的等效节点力，并集成到整体节点力向量中
    [df_number, dummy] = size( gDF ) ;
    for idf = 1:1:df_number
        enf = EquivalentDistPressure( gDF(idf,1:3), gDF(idf,4:6) ) ;
        i = gDF(idf, 1) ;
        j = gDF(idf, 2) ;
        m = gDF(idf, 3) ;
        f( (i-1)*2+1:(i-1)*2+2, 1 ) = f( (i-1)*2+1:(i-1)*2+2, 1 ) + enf( 1:2 ) ;
        f( (j-1)*2+1:(j-1)*2+2, 1 ) = f( (j-1)*2+1:(j-1)*2+2, 1 ) + enf( 3:4 ) ;
        f( (m-1)*2+1:(m-1)*2+2, 1 ) = f( (m-1)*2+1:(m-1)*2+2, 1 ) + enf( 5:6 ) ;
    end
  
    % step5. 处理斜约束边界条件
    [bc3_number,dummy] = size(gBC3) ;
    for ibc=1:1:bc3_number
        n = gBC3( ibc, 1 ) ;
        theta = gBC3( ibc, 4 ) ;
        c = cos(theta) ;
        s = sin(theta) ;
        T = [ c s
             -s c ] ;
        gK((n-1)*2+1:(n-1)*2+2,:) = T*gK((n-1)*2+1:(n-1)*2+2,:) ;
        gK(:,(n-1)*2+1:(n-1)*2+2) = gK(:,(n-1)*2+1:(n-1)*2+2)*transpose(T) ;
        f((n-1)*2+1:(n-1)*2+2,:) = T*f((n-1)*2+1:(n-1)*2+2,:) ;
        gBC1 = [gBC1; gBC3(ibc,1:3)] ;
    end
    
    % step6. 处理第一类约束条件，修改刚度矩阵和节点力向量。采用乘大数法
    [bc1_number,dummy] = size( gBC1 ) ;
    for ibc=1:1:bc1_number
        n = gBC1(ibc, 1 ) ;
        d = gBC1(ibc, 2 ) ;
        m = (n-1)*2 + d ;
        f(m,:) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
        gK(m,m) = gK(m,m) * 1e15 ;
    end
    
    % step7. 求解方程组，得到节点位移向量
    gDelta1 = gK \ f ;
    
    % step8. 把有斜约束的结点位移转换到整体坐标
    gDelta = gDelta1 ;
    for ibc=1:bc3_number
        n = gBC3( ibc, 1 ) ;
        theta = gBC3( ibc, 4 ) ;
        c = cos(theta) ;
        s = sin(theta) ;
        T = [ c s
             -s c ] ;
        gDelta( (n-1)*2+1:n*2,:) = transpose(T)*gDelta( (n-1)*2+1:n*2,: ) ;
    end
    
    % step9. 计算每个单元的结点应力
    gElementStress = zeros( element_number, 8, 3, 2 ) ;
    delta = zeros( 16, 2 ) ;
    for ie = 1:element_number
        xi  = [ -1   1   1  -1   0   1   0  -1 ] ;
        eta = [ -1  -1   1   1  -1   0   1   0 ] ;
        for n=1:8
            B = MatrixB( ie, xi(n), eta(n) ) ;
            D = MatrixD( ie, 1 ) ;
            delta(1:2:15, :) = gDelta( gElement(ie,1:8)*2-1, : ) ;
            delta(2:2:16, :) = gDelta( gElement(ie,1:8)*2, : ) ;
            sigma = D*B*delta ;
            gElementStress( ie, n, :, : ) = sigma ;
        end
    end
return

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

    k = zeros( 16, 16 ) ;
    D = MatrixD( ie, 1 ) ;
    % 3 x 3  高斯积分点和权系数
    %x = [-0.774596669241483,              0.0,0.774596669241483] ;   
    %w = [ 0.555555555555556,0.888888888888889,0.555555555555556] ;  
    % 2 x 2  高斯积分点和权系数
    x = [-0.577350269189626, 0.577350269189626] ;                    
    w = [1, 1] ;                                                     
    for i=1:1:length(x)
        for j=1:1:length(x)
            B = MatrixB( ie, x(i), x(j) ) ;
            J = Jacobi( ie, x(i), x(j) ) ;
            k = k + w(i)*w(j)*transpose(B)*D*B*det(J) ;   
        end
    end
return

function D = MatrixD( ie, opt )
%  计算单元的弹性矩阵D
%  输入参数：
%     ie --------- 单元号
%    opt --------- 问题选项
%                   1 ==>  平面应力
%                   2 ==>  平面应变
%  返回值：
%     D  --------- 弹性矩阵D

    global gElement gMaterial 
    E  = gMaterial( gElement(ie, 9), 1 ) ;   % 弹性模量
    mu = gMaterial( gElement(ie, 9), 2 ) ;   % 泊松比 
    if opt == 1   % 平面应力的弹性常数
        A1 = mu ;                                
        A2 = (1-mu)/2 ;                          
        A3  = E/(1-mu^2) ;                      
    else          % 平面应变的弹性常数
       A1 = mu/(1-mu) ;                        
        A2 = (1-2*mu)/2/(1-mu) ;                
        A3 = E*(1-mu)/(1+mu)/(1-2*mu) ;         
    end
    D = A3* [  1  A1   0
              A1   1   0
               0   0  A2] ;
return

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

    [N_x,N_y] = N_xy( ie, xi, eta );