《有限元基础编程百科全书》| 平面Timoshenko梁单元

上一节中讨论了只承受轴力的杆单元，我们现在来研究一下不仅能承受轴力，还能承受剪力、弯矩的梁单元。同样都为杆系单元，只因承受的的荷载类型不一样，名称叫法也随之不同，理论内容也体现出较多不同。本次我们讨论单元类型有：

• 平面Euler-Bernoulli梁单元
• 空间Euler-Bernoulli梁单元
• 平面Timoshenko梁单元

空间的Timoshenko梁单元程序精度有点问题，就不在这里展示啦。有能力有兴趣的小伙伴可自行编制。本节的参考书籍曾攀^[1]徐荣桥^[2]崔济东^[3]Ferreira^[4]Logan^[5]均已列出，想要深入研究的童鞋，可翻阅相关书籍文献进行深耕。

以上内容分三次推文进行分享，本期分享的是平面Timoshenko梁单元

主要内容有：

• 单元描述
• “叠加原则”计算单元刚度矩阵
• 参照Abaqus帮助文档对剪切部分刚度矩阵修正
• 相关代码
• 数值案例
• 如何定义单元集？

本期完整版代码已存放至《有限元基础编程百科全书》

Timoshenko梁在受力的过程中考虑了剪切变形和转动惯量的影响，即：变形后横截面的平面不在与轴线保持垂直，见下图。

 

Abaqus中的单元号：B21、B22、B31、B31OS、B32、B32OS、PIPE21、PIPE22、PIPE31、PIPE32均是Timoshenko梁。

单元描述

单元自由度表示和Euler-Bernoulli梁是一样的，重点在于刚度矩阵的求解上。

单元刚度矩阵计算

Timoshenko梁单元的刚度矩阵也是通过叠加原则进行推导，只是外加考虑剪切刚度矩阵，除此之外，我们为了对标Abaqus，可以仿照Abaqus对剪切部分的刚度矩阵进行修正，相关理论可点击查看：

 

具体啥意思呢？我来自己总结下，不一定对，仅供参考哈！

其中，    表示修正后的剪切刚度矩阵，    表示未被修正的剪切刚度矩阵，    表示的是一个无量纲因子，预防梁的剪切刚度过大。

其中，    表示截面面积，对于一阶单元    为1.0，对于二阶单元    为    ，SCF默认值为0.25，    为截面惯性矩，    为梁长。

1. 对应于上图中的节点位移         
   轴向位移，直接应用杆单元的刚度矩阵，见式：    
2. 对应于上图中的节点位移         
   该情况是梁纯弯曲情况，其刚度矩阵为：    

   其中，      为截面惯性矩。
3. 对应于上图中的节点位移         
   该情况是梁考虑剪切变形，刚度矩阵可表示为：    

   其中，      为剪切模量，      为横截面面积，      为梁长，      为剪应力不均匀系数，对于      的取值可参考Abaqus帮助文档给出的表格，其中矩形截面，      取0.85。

 

最终叠加后的单元刚度矩阵为：

相关代码

function Ke = Timoshenko2DElementKe(prop,R,L,icoord)
% 单元刚度矩阵
   E = prop(1);v = prop(2);A = prop(3); I = prop(4);
   G = E*0.5/(1+v);
   k = 0.85;% shear factor
   ka3 = k*G*A/L;
   SCF = 0.25;
   xi = 1.0; % 1.0 for first-order elements and 10−4 for second-order elements
   fpa = 1/(1+xi*SCF*L^2*A/12/I);
   KA3 = fpa*ka3;
   KA1 = E*A/L;
   KA2 = E*I/L;
   % Axis
   kea = KA1*[1 0 0 -1 0 0
              0 0 0 0 0 0
              0 0 0 0 0 0
              -1 0 0 1 0 0
              0 0 0 0 0 0
              0 0 0 0 0 0];
   % Bend
   keb = KA2*[0 0 0 0 0 0
               0 0 0 0 0 0 
               0 0 1 0 0 -1
               0 0 0 0 0 0
               0 0 0 0 0 0 
               0 0 -1 0 0 1];
   % Shear(reduced integration)
   kes = KA3*[0 0 0 0 0 0
              0 1 L/2 0 -1 L/2
              0 L/2 L^2/4 0 -L^2/2 L^2/4
              0 0 0 0 0 0
              0 -1 -L/2 0 1 -L/2
              0 L/2 L^2/4 0 -L/2 L^2/4];
   ke = kea + kes + keb;
    switch icoord
        case 1
            Ke = R'*ke*R;
        case 2
            Ke = ke;
    end
   
end

数值案例

本次的案例将采用Abaqus对于Timoshenko梁剪切修正的方式进行数值编程，将单元刚度矩阵与Abaqus导出的单元刚度矩阵进行对标，位移场结果进行对标。

 

 

以1号单元刚度矩阵为例：

 

 

位移场结果对比如下，与Abaqus精度几乎吻合！

自编程序：

     

<<< 左右滑动见更多 >>>

Abaqus：

     

<<< 左右滑动见更多 >>>

如何定义单元集？

本次所用的模型用到了两种材料，这个如何在程序中进行开展呢？首先，我修改了关键词函数，使之可以读取单元集 合信息，关键词定义为：

*Eset1
6,  7,  8, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 31,32, 33
*Eset2
1,  2,  3,  4,  5,  9, 10, 11, 12, 13, 26, 27, 28, 29, 30
*Material
3.0e10,   0.2, 0.16,0.0021333
3.0e10,   0.2, 0.08,0.0010667

讲单元集信息存储为cell类型的变量：

function [Node, Element, Eset, Material, Load, Constr] = MFEATimoshenko2Ddata(filename)
    addpath("input\")
    fid = fopen(strcat(filename,'.txt'), 'r');
    Node = [];
    Element = [];
    Material = [];
    Load = [];
    Constr = [];
    Eset = {}; % Initialize Eset cell array
    
    isNodeSection = false;
    isElementSection = false;
    isMaterialSection = false;
    isLoadSection = false;
    isConstrSection = false;
    isEsetSection = false;
    
    currentEsetIndex = 0; % Initialize current Eset index

    while ~feof(fid)
        line = fgetl(fid);
    
        if isempty(line)
            continue;
        end
    
        % Section checks
        if strncmp(line, '*Node', 5)
            updateSectionFlags('node');
        elseif strncmp(line, '*Element', 8)
            updateSectionFlags('element');
        elseif strncmp(line, '*Material', 9)
            updateSectionFlags('material');
        elseif strncmp(line, '*Load', 5)
            updateSectionFlags('load');
        elseif strncmp(line, '*Constr', 7)
            updateSectionFlags('constr');
        elseif strncmp(line, '*Eset', 5)
            updateSectionFlags('eset');
            currentEsetIndex = currentEsetIndex + 1; % Increment Eset index for a new set
            Eset{currentEsetIndex} = []; % Initialize new row for Eset
            continue;
        end
    
        % Data parsing
        parseData(line);
    end
    
    fclose(fid);
end

在计算单元刚度矩阵时，对于不同的单元集赋予不同的材料属性：

% 为梁单元和柱单元分别定义材料属性
BeamMaterial = Material(2,:);
PillarMaterial = Material(1,:);

% 根据 Eset{1} 和 Eset{2} 获取梁单元和柱单元的索引
BeamElementIndices = Eset{2};
PillarElementIndices = Eset{1};
for iEle = 1:EleCount
    % 检查当前单元属于哪种类型，并选择相应的材料属性
    if ismember(iEle, BeamElementIndices)
        CurrentMaterial = BeamMaterial;
    elseif ismember(iEle, PillarElementIndices)
        CurrentMaterial = PillarMaterial;
    else
        error('未知单元类型');
    end
    % 获取当前单元的节点索引
    n1 = Element(iEle, 1); n2 = Element(iEle, 2);
    R = CoordTransform([x(n1), x(n2)], [y(n1), y(n2)], BarLength(iEle));
    ElementStiffnessMatrix = Timoshenko2DElementKe(CurrentMaterial, R, BarLength(iEle), 1);
    % 确定节点1和节点2的自由度范围
    n1_dofs = (n1 - 1) * Dof + (1:Dof);
    n2_dofs = (n2 - 1) * Dof + (1:Dof);
    % 将节点1和节点2的自由度合并为一个向量
    ElementNodeDOF = [n1_dofs, n2_dofs];
    KKG(ElementNodeDOF, ElementNodeDOF) = KKG(ElementNodeDOF, ElementNodeDOF) + ElementStiffnessMatrix;
end

在程序中可以导出每个单元的刚度矩阵，用于和Abaqus对比，采用以下程序段：

matObj = matfile('ElementStiffnessMatrix.mat', 'Writable', true);
varName = sprintf('ElementStiffnessMatrix%d', iEle); % 动态创建变量名
matObj.(varName) = ElementStiffnessMatrix;

本章小结

经过多天的打磨，终将杆系结构章节更新完毕，该章节采用全新的彩色插图，努力使读者眼前一亮，每个单元之前均给出了整套代码和视频讲解链接，希望可以带着小白们步入有限元数值世界的大门。

杆系单元的刚度矩阵均给出了具体形式，用户在编程时可以直接赋予具体的数值，相对较为简单，刚度矩阵的组成也严格遵循"叠加原则"，相对较为容易理解，紧紧把握住"坐标变换"和"叠加原则"这两个重要概念。

在下一个章节中，我将分享平面问题的有限元问题，届时会带来等参思想、各种应力、应变概念、数值积分方案求解刚度矩阵、低阶&高阶单元、应力平均、节点场插值等等精彩内容，欢迎读者继续关注！

参考资料