首页/文章/ 详情

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

6月前浏览785


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

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

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

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

主要内容有:

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

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


单元描述

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

 

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

 
 

相应的刚度方程为:

 

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

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

 

看似很麻烦,实则一点也不简单(开玩笑啦哈哈哈哈),可以根据刚度矩阵叠加原则,一点一点叠加上去!

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

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

 

空间坐标变换

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

数值案例

建立下图所示的空间Euler-Bernoulli梁有限元模型,单元节点编码、荷载形式、力学参数均已标注,试求各节点自由度方向上的位移量,单元节点力留作为读者兴趣研究。

 
 
 
 
 

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