有限元基础编程——CST单元

之前的推文中讲述的有限元编程过程中的单元刚度组集过程中，往往是直接照搬现成的单元刚度矩阵，这个做法对于简单的梁单元和杠单元来说是适用的，对于连续体结构有限元分析时，使用直接刚度法显得很“笨重”，单元刚度矩阵往往需要结合最小势能原理推导的公式进行计算，本次分享的时结构平面分析最简单的“入门”单元——CST单元(常应变三角形单元)，延续前几篇的风格，首先介绍基本理论，然后基于Matlab对其进行编程。

理论部分

位移场

CST单元共有三个节点，每个节点两个自由度，单元内任意点的位移可通过形函数    进行插值得到： 

矩阵可表示为：

上图利用等参变换思想将任意三角形规则化，单元形函数    为

雅可比矩阵

几何方程与应变矩阵

    令    ，其中

我们可以从以上式子中看出，应变矩阵    是一个常数矩阵，这也就是其称之为“常应变矩阵”的原因，对于线弹性范围内，其应力矩阵也将是常应力。

单元刚度矩阵

利用最小势能原理可得单元刚度矩阵     

由于B矩阵为常数矩阵，上式可写为：

A为三角形单元的面积，等于    。

单元边界处理

由于计算模型只有两个单元，故进行”降阶“处理，即划行划列，至于单元数量较多的情况用到的置1法与乘大数法将会在以后专门推出一期进行讲解。

理论部分终于写完了，呼呼呼~接下来是代码部分咯，可直接Copy移植哦

模型描述

如上图所示的一个薄板，在右端受集中力$F$作用，其中参数为：    

代码区

function [J B Ae]=CST_B(xy)
% 用于计算B矩阵
    NN=[1 0 -1
        0 1 -1];
    J = NN*xy;
    Ae = 0.5*det(J);
    A = 1/det(J)*[J(2,2)  -J(1,2)  0   0
                    0      0   -J(2,1)  J(1,1)
                -J(2,1)  J(1,1)    J(2,2) -J(1,2)];
    G = [1 0 0 0 -1 0
         0 0 1 0 -1 0
         0 1 0 0 0 -1
         0 0 0 1 0 -1];
    B = A*G;
end

function k = CST_Stiffness(E,NU,t,xe,ye,ID)
% 计算单元刚度矩阵
    if ID == 1
         D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
    elseif ID == 2
         D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
    end

    xy = zeros(3,2);
    for i = 1:3
        xy(i,1)=xe(i);
        xy(i,2)=ye(i);
    end

    k = zeros(6,6);
    [J B Ae]=CST_B(xy);
    k = t*Ae*B'*D*B;
end

function z = Triangle2D3Node_Assembly(KK,k,i,j,m)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号I、j、m
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
DOF(5)=2*m-1;
DOF(6)=2*m;
for n1=1:6
   for n2=1:6
      KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
   end
end
z=KK;

function stress=CST_Stress(E,NU,xe,ye,ele_u,ID)
%该函数计算单元的应力
%输出单元的应力stress(3X1)，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sz
if ID == 1 
   D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
elseif ID == 2
   D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
end
xy = zeros(3,2);
for i = 1:3
    xy(i,1)=xe(i);
    xy(i,2)=ye(i);
end
[J B Ae]=CST_B(xy);
stress = D*B*ele_u;