美文网首页
有限元基础编程——杆单元(附Matlab源码)

有限元基础编程——杆单元(附Matlab源码)

作者: 易木木响叮当 | 来源:发表于2022-07-28 17:26 被阅读0次

    本篇内容转载本人公众号:易木木响叮当,涉及代码可以在后台回复:杆_Matlab,即可自动获取。

    引言

    ”有限的单元,无限的能力“这句话来自清华大学有限元分析公开课曾攀老师的开课语。想要学好有限元这门课,不光要理解理论公式的由来及简单手酸,更要结合实际应用。本栏目将带着大家Step-By-Step基于Matlab语言实现有限元的基础操作,课程代码来自《有限元分析基础教程》——曾攀,并附赠ANSYS命令流文件进行验证Matlab代码正确性。

    有限元“流水线套路”

    • 求解单元刚度

    • 组装整体刚度

    • 未知位移求解

      • 本质是线性方程组求解,求解方法有很多,基于Fortran编写的可以采用JCG开源程序包,基于Matlab编写的可以采用\,默认高斯消去法,也可以使用PCG(预处理共轭梯度法)等等,总之求解线性方程组的方法很多,我们初学者可以先使用最简单的\,进行求解。

      • 对于位移边界的处理,有限元有很多方法:直接消去法置1法拉格朗日乘子法罚函数法等,对于初学者可以先概念最简单的直接消去法入手,等熟悉了有限元基本过程再使用更加进阶的操作。

    • 节点力、应力、应变等求解

    1D杆单元

    这大概是最简单的有限元分析吧,简直每本有限元教材里面都会将之作为入门案例,操作虽然简单,但也是包含了有限元分析的基础步骤。

    案例

    例题.jpg

    函数

    function k=Bar1D2Node_Stiffness(E,A,L)
    k=[E*A/L -E*A/L; -E*A/L E*A/L];
    function z=Bar1D2Node_Assembly(KK,k,i,j)
    DOF(1)=i;
    DOF(2)=j;
    for n1=1:2
      for n2=1:2
          KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
      end
    end
    z=KK;
    

    主程序

    format long
    % 典型例题[2.3(1)]P17
    E1 = 2*10^5;E2 = E1;E3 = E1;
    A3 = 0.06;A2 = 0.5*A3;A1 = A3/3;
    L1 = 0.1;L2 = L1;L3 =L1;
    k1 = Bar1D2Node_Stiffness(E1,A1,L1);
    k2 = Bar1D2Node_Stiffness(E2,A2,L2);
    k3 = Bar1D2Node_Stiffness(E3,A3,L3);
    KK = zeros(4,4);
    KK = Bar1D2Node_Assembly(KK,k1,1,2);
    KK = Bar1D2Node_Assembly(KK,k2,2,3);
    KK = Bar1D2Node_Assembly(KK,k3,3,4);
    % 直接法求解整体刚度矩阵
    k = KK([1:3],[1:3]);
    p = [-100;0;50];
    % 高斯消去法求解线性方程组
    u = k\p
    

    结果所得位移与手算结果一致,可自行验证。

    2D杆单元

    坐标变换

    2D杆单元在编写的时候涉及到由局部坐标系向整体坐标系变换的过程。坐标转换矩阵T为:

    T=\left[ \begin{matrix} \cos \alpha& \sin \alpha& 0& 0\\ 0& 0& \cos \alpha& \sin \alpha\\ \end{matrix} \right]

    刚度矩阵、节点位移由局部坐标系Kq转换到整体坐标系\bar{K}\bar{q}

    \bar{K}=T^TKT \\ \bar{q}=T^Tq

    \bar{K}=\frac{EA}{L}\left[ \begin{matrix} \cos ^2\alpha& \cos \alpha \sin \alpha& -\cos ^2\alpha& -\cos \alpha \sin \alpha\\ \cos \alpha \sin \alpha& \sin ^2\alpha& -\cos \alpha \sin \alpha& -\sin ^2\alpha\\ -\cos ^2\alpha& -\cos \alpha \sin \alpha& \cos ^2\alpha& \cos \alpha \sin \alpha\\ -\cos \alpha \sin \alpha& -\sin ^2\alpha& \cos \alpha \sin \alpha& \sin ^2\alpha\\ \end{matrix} \right]

    进行应力、节点力计算时,位移也应该由局部坐标系转换到整体坐标系中,由弹性力学中的物理方程,有1D问题的应力:

    \sigma =E\varepsilon =EBq=\frac{E}{L}\left[ \begin{matrix} -1& 1\\ \end{matrix} \right] Tq

    F=\sigma A

    案例

    杆单元案例

    函数

    function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)
    %该函数计算单元的刚度矩阵
    %输入弹性模量E,横截面积A
    %输入第一个节点坐标(x1,y1),第二个节点坐标(x2,y2),角度alpha(单位是度)
    %输出单元刚度矩阵k(4X4)。
    
    L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
    x=alpha*pi/180;
    C=cos(x);
    S=sin(x);
    k=E*A/L*[C*C C*S -C*C -C*S
    C*S S*S -C*S -S*S
    -C*C -C*S C*C C*S
    -C*S -S*S C*S S*S];
    function z = Bar2D2Node_Assembly(KK,k,i,j)
    %该函数进行单元刚度矩阵的组装
    %输入单元刚度矩阵k,单元的节点编号i、j
    %输出整体刚度矩阵KK
    
    DOF(1)=2*i-1;
    DOF(2)=2*i;
    DOF(3)=2*j-1;
    DOF(4)=2*j;
    for n1=1:4
      for n2=1:4
          KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
      end
    end
    z=KK;
    

    主程序

    E=2.95e11;A=0.0001;
    x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;
    alpha1=0;alpha2=90;alpha3=atan(0.75)*180/pi;
    k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2,alpha1)
    k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3,alpha2) 
    k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3)
    k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3,alpha1) 
    %建立整体刚度方程
    %由于该结构共有4个节点,因此,结构总的刚度矩阵为KK(8×8),先对K清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。
    KK=zeros(8,8);
    KK=Bar2D2Node_Assembly (KK,k1,1,2);
    KK=Bar2D2Node_Assembly (KK,k2,2,3);
    KK=Bar2D2Node_Assembly (KK,k3,1,3);
    KK=Bar2D2Node_Assembly (KK,k4,4,3)
    %边界条件的处理及刚度方程求解
    k=KK([3,5,6],[3,5,6])
    p=[20000;0;-25000]
    u=k\p
    

    ANSYS命令流

    / PREP7              !进入前处理 
    /PLOPTS,DATE,0      !设置不显示日期和时间
    !=====设置单元、材料,生成节点及单元
    ET,1,LINK1            !选择单元类型
    UIMP,1,EX, , ,2.95e11,    !给出材料的弹性模量
    R,1,1e-4,              !给出实常数(横截面积)
    N,1,0,0,0,              !生成1号节点,坐标(0,0,0) 
    N,2,0.4,0,0,            !生成2号节点,坐标(0.4,0,0)
    N,3,0.4,0.3,0,            !生成3号节点,坐标(0.4,0.3,0)
    N,4,0,0.3,0,            !生成4号节点,坐标(0,0.3,0)
    E,1,2                  !生成1号单元(连接1号节点和2号节点) 
    E,2,3                  !生成2号单元(连接2号节点和3号节点)
    E,1,3                  !生成3号单元(连接1号节点和3号节点)
    E,4,3                  !生成4号单元(连接4号节点和3号节点)
    FINISH                !前处理结束
    !=====在求解模块中,施加位移约束、外力,进行求解
    /SOLU                  !进入求解状态(在该状态可以施加约束及外力)
    D,1,ALL                !将1号节点的位移全部固定
    D,2,UY,                !将2号节点的Y方向位移固定
    D,4,ALL                !将4号节点的位移全部固定
    F,2,FX,20000,            !在2号节点处施加X方向的力(20000)
    F,3,FY,-25000,            !在3号节点处施加Y方向的力(-25000)
    SOLVE                  !进行求解
    FINISH                  !结束求解状态
    !=====进入一般的后处理模块
    /POST1                  !进入后处理
    PLDISP,1                  !显示变形状况
    FINISH                  !结束后处理
    
    

    相关文章

      网友评论

          本文标题:有限元基础编程——杆单元(附Matlab源码)

          本文链接:https://www.haomeiwen.com/subject/eaeqwrtx.html