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

6月前浏览1380

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

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

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

以上内容分三次推文进行分享，本期分享的是2D Euler-Bernoulli梁单元

主要内容有：

单元描述
单元刚度矩阵计算
等效节点荷载
相关代码
数值案例

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

Euler-Bernoulli梁不考虑剪切变形和转动惯量的影响，基于两个假设：

变形前垂直梁中心线的横截面，变形后仍为平面（刚性横截面假定）；
变形后横截面的平面仍与变形后的轴线垂直，详见下图

Euler-Bernoulli梁适合跨高比较大的梁结构，在总变形部分中，剪切变形占比重较小，一般可以忽略不计。在Abaqus中B23、B23H、B33、和B33H单元号表示Euler-Bernoulli梁单元。

单元描述

局部坐标系下的梁单元每个节点具有6个自由度，见下图

上图所示的平面梁单元节点位移列向量和节点力列向量分别为:

相应的刚度方程为：

单元刚度矩阵计算

若不计轴向位移的自由度，单元刚度矩阵可表示为：

考虑轴向位移后，单元刚度矩阵可根据叠加原则，将局部坐标系下的杆单元刚度矩阵，根据自由度号添加至上式，形成局部坐标系下的平面梁单元单元刚度矩阵：

等效节点荷载

有时候的荷载形式并不是直接施加在节点上，而是施加在单元面上，无法直接参加节点外荷载向量组装，这时候该怎么办呢？可以通过非节点荷载虚功相等的原则将荷载等效到节点上，即：

下图形象地表达了承受均布荷载平面梁单元如何转化为节点荷载:

相关代码

function Ke = Bernoulli2DElementKe(prop,R,L,icoord)
% 单元刚度矩阵
   E = prop(1);I = prop(2); A = prop(3);
   ke=[E*A/L       0             0          -E*A/L       0                0
       0        12*E*I/L^3      6*E*I/L^2     0       -12*E*I/L^3   6*E*I/L^2;
       0        6*E*I/L^2       4*E*I/L       0        -6*E*I/L^2 2*E*I/L ;
       -E*A/L       0             0          E*A/L       0                0
       0        -12*E*I/L^3     -6*E*I/L^2     0       12*E*I/L^3   -6*E*I/L^2;
       0        6*E*I/L^2       2*E*I/L       0        -6*E*I/L^2 4*E*I/L ;];
   
    switch icoord
        case 1
            Ke = R'*ke*R;
        case 2
            Ke = ke;
    end
   
end

function enf = EquivalentNodeForce(Node,Element,iEle,p1, p2, idof )
%   计算线性分布荷载的等效节点力
%   输入参数:
%      ie  -----  单元号
%      p1  -----  第一个节点上的分布力集度值
%      p2  -----  第二个节点上的分布力集度值
%      idof  ---  分布力的种类，它可以是下面几种
%                  1 ---  分布轴向力
%                  2 ---  分布横向力
%                  3 ---  分布弯矩
%   返回值:
%      enf -----  整体坐标系下等效节点力向量
    enf = zeros( 6, 1 ) ;                       % 定义 6x1 的等效节点力向量
    x = Node(:,1);y = Node(:,2);
    BarLength=BarsLength(x,y,Element);
    L = BarLength(iEle);
    n1=Element(iEle,1);n2=Element(iEle,2);
    R = CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],L);
    switch idof 
    case 1     %  分布轴向力 
        enf( 1 ) = (2*p1+p2)*L/6 ;
        enf( 4 ) = (p1+2*p2)*L/6 ;
    case 2     %  分布横向力
        enf( 2 ) = (7*p1+3*p2)*L/20 ;
        enf( 3 ) = (3*p1+2*p2)*L^2/60 ;
        enf( 5 ) = (3*p1+7*p2)*L/20 ;
        enf( 6 ) = -(2*p1+3*p2)*L^2/60 ;
    case 3     %  分布弯矩
        enf( 2 ) = -(p1+p2)/2 ;
        enf( 3 ) = (p1-p2)*L/12 ;
        enf( 5 ) = (p1+p2)/2 ;
        enf( 6 ) = -(p1-p2)*L/12 ;
    end
    enf = R' * enf  ;                         %  把等效节点力转换到整体坐标下
end