《有限元基础编程百科全书》| 平面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， 为截面惯性矩， 为梁长。
对应于上图中的节点位移
轴向位移，直接应用杆单元的刚度矩阵，见式：对应于上图中的节点位移
该情况是梁纯弯曲情况，其刚度矩阵为： 其中， 为截面惯性矩。对应于上图中的节点位移
该情况是梁考虑剪切变形，刚度矩阵可表示为： 其中， 为剪切模量， 为横截面面积， 为梁长， 为剪应力不均匀系数，对于 的取值可参考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;
本章小结
经过多天的打磨，终将杆系结构章节更新完毕，该章节采用全新的彩色插图，努力使读者眼前一亮，每个单元之前均给出了整套代码和视频讲解链接，希望可以带着小白们步入有限元数值世界的大门。
杆系单元的刚度矩阵均给出了具体形式，用户在编程时可以直接赋予具体的数值，相对较为简单，刚度矩阵的组成也严格遵循"叠加原则"，相对较为容易理解，紧紧把握住"坐标变换"和"叠加原则"这两个重要概念。
在下一个章节中，我将分享平面问题的有限元问题，届时会带来等参思想、各种应力、应变概念、数值积分方案求解刚度矩阵、低阶&高阶单元、应力平均、节点场插值等等精彩内容，欢迎读者继续关注！
参考资料
曾攀: 《有限元分析基础教程》
[2]徐荣桥: 《结构分析的有限元法与MATLAB程序设计》
[3]崔济东: 《有限单元法——编程与软件应用》
[4]Ferreira: 《MATLAB-codes-for-finite-element-analysis》
[5]Logan: 《A-first-course-in-the-finite-element-method》
