题外话:因为微 信的推荐机制变动,有可能大家不会第一时间看到我的文章,请大家给我的公 众号标上⭐,以免错过好资源。
上一节中讨论了只承受轴力的杆单元,我们现在来研究一下不仅能承受轴力,还能承受剪力、弯矩的梁单元。同样都为杆系单元,只因承受的的荷载类型不一样,名称叫法也随之不同,理论内容也体现出较多不同。本次我们讨论单元类型有:
空间的Timoshenko梁单元程序精度有点问题,就不在这里展示啦。有能力有兴趣的小伙伴可自行编制。本节的参考书籍曾攀[1]徐荣桥[2]崔济东[3]Ferreira[4]Logan[5]均已列出,想要深入研究的童鞋,可翻阅相关书籍文献进行深耕。
以上内容分三次推文进行分享,本期分享的是平面Timoshenko梁单元
主要内容有:
本期完整版代码已存放至《有限元基础编程百科全书》
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
具体啥意思呢?我来自己总结下,不一定对,仅供参考哈!
其中,
其中,
最终叠加后的单元刚度矩阵为:
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》