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

6月前浏览775

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

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

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

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

主要内容有：

单元描述
“叠加原则”计算单元刚度矩阵
空间坐标转换
相关代码
数值案例

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

单元描述

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

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

相应的刚度方程为：

“叠加原则”计算单元刚度矩阵

空间梁单元刚度矩阵为一个的矩阵，见下图

看似很麻烦，实则一点也不简单（开玩笑啦哈哈哈哈），可以根据刚度矩阵叠加原则，一点一点叠加上去！

对应于上图中的节点位移
轴向位移，直接应用杆单元的刚度矩阵：
对应于上图中的节点位移
扭转角，其刚度矩阵为：
其中，为横截面的扭转惯性矩，为剪切模量。
对应于上图中的节点位移
该情况是梁在平面内纯弯曲情况，对应于：
其中为绕着轴的中性轴的惯性矩。
对应于上图中的节点位移该情况是梁在平面内纯弯曲情况：
其中为绕着轴的中性轴的惯性矩。

最终组装的单元刚度矩阵为下面这个样子：

空间坐标变换

3D的坐标变换比较复杂，理论拎不清楚（能力确实有限~），我就直接放代码了！

相关代码

function R = CoordTransform(x,y,z,L)
% 坐标转换
    x1 = x(1);x2 = x(2);y1 = y(1);y2 = y(2);z1 = z(1);z2 = z(2);
    if x1 == x2 && y1 == y2
        if z2 > z1
            Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
        else
            Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
        end
    else
        CXx = (x2-x1)/L;
        CYx = (y2-y1)/L;
        CZx = (z2-z1)/L;
        D = sqrt(CXx*CXx + CYx*CYx);
        CXy = -CYx/D;
        CYy = CXx/D;
        CZy = 0;
        CXz = -CXx*CZx/D;
        CYz = -CYx*CZx/D;
        CZz = D;
        Lambda = [CXx CYx CZx ;CXy CYy CZy ;CXz CYz CZz];
    end
    R = [Lambda zeros(3,9); zeros(3) Lambda zeros(3,6);
    zeros(3,6) Lambda zeros(3);zeros(3,9) Lambda];
end

function Ke = Bernoulli3DElementKe(prop,R,L,icoord)
% 单元刚度矩阵
    E = prop(1);A = prop(2);Iy = prop(3);Iz = prop(4);G = prop(5);J = prop(6); 
    k1 = E*A/L;
    k2 = 12*E*Iz/(L*L*L);
    k3 = 6*E*Iz/(L*L);
    k4 = 4*E*Iz/L;
    k5 = 2*E*Iz/L;
    k6 = 12*E*Iy/(L*L*L);
    k7 = 6*E*Iy/(L*L);
    k8 = 4*E*Iy/L;
    k9 = 2*E*Iy/L;
    k10 = G*J/L;
    ke = [k1 0 0 0 0 0 -k1 0 0 0 0 0;
          0 k2 0 0 0 k3 0 -k2 0 0 0 k3;
          0 0 k6 0 -k7 0 0 0 -k6 0 -k7 0;
          0 0 0 k10 0 0 0 0 0 -k10 0 0;
          0 0 -k7 0 k8 0 0 0 k7 0 k9 0;
          0 k3 0 0 0 k4 0 -k3 0 0 0 k5;
         -k1 0 0 0 0 0 k1 0 0 0 0 0;
          0 -k2 0 0 0 -k3 0 k2 0 0 0 -k3;
          0 0 -k6 0 k7 0 0 0 k6 0 k7 0;
          0 0 0 -k10 0 0 0 0 0 k10 0 0;
          0 0 -k7 0 k9 0 0 0 k7 0 k8 0;
          0 k3 0 0 0 k5 0 -k3 0 0 0 k4];
   
    switch icoord
        case 1
            Ke = R'*ke*R;
        case 2
            Ke = ke;
    end
   
end