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

11天前浏览249

题外话：因为微信的推荐机制变动，有可能大家不会第一时间看到我的文章，请大家给我的公众号标上⭐，以免错过好资源。

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

平面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对剪切部分的刚度矩阵进行修正，相关理论可点击查看：Abaqus官方文档：
https://help.3ds.com/2024/English/DSSIMULIA_Established/SIMACAEELMRefMap/simaelm-c-beamelem.htm?contextscope=all&id=a2e4f2c8942342699ebd42cf79509e1c

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

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

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

对应于上图中的节点位移
轴向位移，直接应用杆单元的刚度矩阵，见式：
对应于上图中的节点位移
该情况是梁纯弯曲情况，其刚度矩阵为：
其中，为截面惯性矩。
对应于上图中的节点位移
该情况是梁考虑剪切变形，刚度矩阵可表示为：
其中，为剪切模量，为横截面面积，为梁长，为剪应力不均匀系数，对于的取值可参考Abaqus帮助文档给出的表格，其中矩形截面，取0.85。

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

数值案例

本次的案例将采用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
% 检查当前单元属于哪种类型，并选择相应的材料属性
ifismember(iEle, BeamElementIndices)
        CurrentMaterial = BeamMaterial;
elseifismember(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;