首页/文章/ 详情

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

5月前浏览1328

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

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

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

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

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

主要内容有:

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

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


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

 

Abaqus中的单元号:B21B22B31B31OSB32B32OSPIPE21PIPE22PIPE31PIPE32均是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,    为截面惯性矩,    为梁长。

  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;

本章小结

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

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

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

参考资料

[1]

曾攀: 《有限元分析基础教程》

[2]

徐荣桥: 《结构分析的有限元法与MATLAB程序设计》

[3]

崔济东: 《有限单元法——编程与软件应用》

[4]

Ferreira: 《MATLAB-codes-for-finite-element-analysis》

[5]

Logan: 《A-first-course-in-the-finite-element-method》

以上推文已更新至《有限元基础编程百科全书》PDF文档中,扫描下方二维码或在后台回复“星球”,加入知识星球后,查看置顶文章即可获得《有限元基础编程百科全书》最新更新版本。



来源:易木木响叮当
ACTAbaqusMATLABADS理论材料
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-05-11
最近编辑:5月前
易木木响叮当
硕士 有限元爱好者
获赞 212粉丝 243文章 342课程 2
点赞
收藏
作者推荐
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈