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

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

理论部分

位移场

有限元基础编程——CST单元的图1

CST单元共有三个节点，每个节点两个自由度，单元内任意点的位移可通过形函数 $N\left( \eta ,\xi \right)$ 进行插值得到：

$\begin{cases} u=N_1u_1+N_2u_2+N_3u_3\\ v=N_1v_1+N_2v_2+N_3v_3\\ \end{cases}$

矩阵可表示为：

$\left[ \begin{array}{c} u\\ v\\ \end{array} \right] =\left[ \begin{matrix} N_1& 0& N_2& 0& N_3& 0\\ 0& N_1& 0& N_2& 0& N_3\\ \end{matrix} \right] \left[ \begin{array}{c} u_1\\ v_1\\ u_2\\ v_2\\ u_3\\ v_3\\ \end{array} \right]$

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

$\begin{cases} N_1=\xi\\ N_2=\eta\\ N_3=1-\xi -\eta\\ \end{cases}$

雅可比矩阵

$J=\left[ \begin{matrix} \frac{\partial x}{\partial \xi}& \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta}& \frac{\partial y}{\partial \eta}\\ \end{matrix} \right] =\left[ \begin{matrix} \sum_{i=1}^3{\frac{\partial N_i}{\partial \xi}x_i}& \sum_{i=1}^3{\frac{\partial N_i}{\partial \xi}y_i}\\ \sum_{i=1}^3{\frac{\partial N_i}{\partial \eta}x_i}& \sum_{i=1}^3{\frac{\partial N_i}{\partial \eta}y_i}\\ \end{matrix} \right]$

$J^{-1}=\frac{1}{\det J}\left[ \begin{matrix} J_{22}& -J_{12}\\ -J_{21}& J_{11}\\ \end{matrix} \right]$

几何方程与应变矩阵

$\varepsilon =\left[ B \right] q^e$

$\varepsilon =\left[ \begin{array}{c} \varepsilon _x\\ \varepsilon _y\\ \gamma _{xy}\\ \end{array} \right] =\frac{1}{\det J}\left[ \begin{matrix} J_{22}& -J_{12}& 0& 0\\ 0& 0& -J_{21}& J_{11}\\ -J_{21}& J_{11}& J_{22}& -J_{12}\\ \end{matrix} \right] \left[ \begin{array}{c} \frac{\partial u}{\partial \xi}\\ \frac{\partial u}{\partial \eta}\\ \frac{\partial v}{\partial \xi}\\ \frac{\partial u}{\partial \eta}\\ \end{array} \right]$

$=\frac{1}{\det J}\left[ \begin{matrix} J_{22}& -J_{12}& 0& 0\\ 0& 0& -J_{21}& J_{11}\\ -J_{21}& J_{11}& J_{22}& -J_{12}\\ \end{matrix} \right] \left[ \begin{matrix} \frac{\partial N_1}{\partial \xi}& 0& \frac{\partial N_2}{\partial \xi}& 0& \frac{\partial N_3}{\partial \xi}& 0\\ \frac{\partial N_1}{\partial \eta}& 0& \frac{\partial N_2}{\partial \eta}& 0& \frac{\partial N_3}{\partial \eta}& 0\\ 0& \frac{\partial N_1}{\partial \xi}& 0& \frac{\partial N_1}{\partial \xi}& 0& \frac{\partial N_1}{\partial \xi}\\ 0& \frac{\partial N_2}{\partial \eta}& 0& \frac{\partial N_2}{\partial \eta}& 0& \frac{\partial N_2}{\partial \eta}\\ \end{matrix} \right] \left[ \begin{array}{c} u_1\\ v_1\\ u_2\\ v_2\\ u_3\\ v_3\\ \end{array} \right]$ 令 $\left[ B \right] =\left[ A \right] \left[ G \right]$ ，其中

$\left[ A \right] =\frac{1}{\det J}\left[ \begin{matrix} J_{22}& -J_{12}& 0& 0\\ 0& 0& -J_{21}& J_{11}\\ -J_{21}& J_{11}& J_{22}& -J_{12}\\ \end{matrix} \right]$

$\left[ G \right] =\left[ \begin{matrix} \frac{\partial N_1}{\partial \xi}& 0& \frac{\partial N_2}{\partial \xi}& 0& \frac{\partial N_3}{\partial \xi}& 0\\ \frac{\partial N_1}{\partial \eta}& 0& \frac{\partial N_2}{\partial \eta}& 0& \frac{\partial N_3}{\partial \eta}& 0\\ 0& \frac{\partial N_1}{\partial \xi}& 0& \frac{\partial N_1}{\partial \xi}& 0& \frac{\partial N_1}{\partial \xi}\\ 0& \frac{\partial N_2}{\partial \eta}& 0& \frac{\partial N_2}{\partial \eta}& 0& \frac{\partial N_2}{\partial \eta}\\ \end{matrix} \right]$

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

单元刚度矩阵

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

$k^e=\int_{\varOmega ^e}{B^TDBd\varOmega}=t\int_{A^e}{B^TDB\,\,dA}$

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

$k^e=B^TDB\,\,t\,\,A$

A为三角形单元的面积，等于 $\frac{1}{2}\det J$ 。

单元边界处理

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

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

模型描述

有限元基础编程——CST单元的图2

如上图所示的一个薄板，在右端受集中力$F$作用，其中参数为： $E=1\times 10^7Pa\text{，}\mu =1/3\text{，}t=0.1m\text{，}F=1\times 10^5N$

代码区

% 篇幅原因只展示子函数，主程序可在公众号：
% 易木木响叮当，后台回复“CST单元”，即可自动获取
% --------------------------------------
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;

有限元基础编程：CST单元实现方法详解