[FEM-7]杆单元分析的MATLAB程序

---

1 1D杆单元的有限元分析程序(Bar1D2Node)

最简单的线性杆单元程序应包括单元刚度矩阵、单元组装、单元应力等几个基本计算程序。

1.1 单元刚度矩阵

Element_K_Array = Bar1D2Node_Stiffness(Mat_E, Geo_A, Geo_L)

• ToDo：计算单元的刚度矩阵
• Input：弹性模量Mat_E、横截面积Geo_A、长度Geo_L
• Output：单元刚度矩阵Element_K

function Element_K_Array = Bar1D2Node_Stiffness(Mat_E, Geo_A, Geo_L)
    Element_K_Array = [Mat_E * Geo_A / Geo_L, -Mat_E * Geo_A / Geo_L; ...
     -Mat_E * Geo_A / Geo_L, Mat_E * Geo_A / Geo_L];
end

图1 一维三连杆结构

上图为一维三连杆结构，通过程序得到每个杆件的单元刚度方程分别为

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{(1)}=\left[\begin{array}{cc}{K}_{11}^{(1)}& K_{12}^{(1)}\\ K_{21}^{(1)}& K_{22}^{(1)}\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{(2)}=\left[\begin{array}{cc}{K}_{22}^{(2)}& K_{23}^{(2)}\\ K_{32}^{(2)}& K_{33}^{(2)}\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{(3)}=\left[\begin{array}{cc}{K}_{33}^{(3)}& K_{34}^{(3)}\\ K_{43}^{(3)}& K_{44}^{(3)}\end{array}\right]\end{array}\end{array}$$

1.2 整体刚度矩阵

Total_K_Array = Bar1D2Node_Assembly(Total_K_Array, Element_K_Array, ElementNode_1, ElementNode_2)

• ToDo：计算整体刚度矩阵
• Input：单元刚度矩阵Element_K_Array、单元两端节点编号ElementNode_1、ElementNode_2
• Output：单元刚度矩阵Total_K_Array

function Total_K_Array = Bar1D2Node_Assembly(Total_K_Array, Element_K, ElementNode_1, ElementNode_2)
    Element_DOF_Array(1) = ElementNode_1;
    Element_DOF_Array(2) = ElementNode_2;
    for ii = 1 : 2
        for jj = 1 : 2
            Total_K_Array(Element_DOF_Array(ii), Element_DOF_Array(jj)) = Total_K_Array(Element_DOF_Array(ii), Element_DOF_Array(jj)) + Element_K_Array(ii, jj);
        end
    end
    Total_K_Array = Total_K_Array;
end

1.2.1 组装第1根杆件

通过程序，初始时Total_K_Array为$4\times 4$的空矩阵，当组装第一根杆件时，有

Total_K_Array = Bar1D2Node_Assembly(Total_K_Array, Element_K_Array_1, 1, 2);

此时，

Element_DOF_Array(1) = 1;
Element_DOF_Array(2) = 2;

循环有，

ii = 1, jj = 1
Total_K_Array(1, 1) = Total_K_Array(1, 1) + Element_K_Array_1(1, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 1, jj = 2
Total_K_Array(1, 2) = Total_K_Array(1, 2) + Element_K_Array_1(1, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 2, jj = 1
Total_K_Array(2, 1) = Total_K_Array(2, 1) + Element_K_Array_1(2, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 2, jj = 2
Total_K_Array(2, 2) = Total_K_Array(2, 2) + Element_K_Array_1(2, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

1.2.2 组装第2根杆件

Total_K_Array = Bar1D2Node_Assembly(Total_K_Array, Element_K_Array_2, 2, 3);

此时，

Element_DOF_Array(1) = 2;
Element_DOF_Array(2) = 3;

循环有，

ii = 1, jj = 1
Total_K_Array(2, 2) = Total_K_Array(2, 2) + Element_K_Array_2(1, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 1, jj = 2
Total_K_Array(2, 3) = Total_K_Array(2, 3) + Element_K_Array_2(1, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 2, jj = 1
Total_K_Array(3, 2) = Total_K_Array(3, 2) + Element_K_Array_2(2, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 2, jj = 2
Total_K_Array(3, 3) = Total_K_Array(3, 3) + Element_K_Array_2(2, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

1.2.3 组装第3根杆件

Total_K_Array = Bar1D2Node_Assembly(Total_K_Array, Element_K_Array_3, 3, 4);

此时，

Element_DOF_Array(1) = 3;
Element_DOF_Array(2) = 4;

循环有，

ii = 1, jj = 1
Total_K_Array(3, 3) = Total_K_Array(3, 3) + Element_K_Array_3(1, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}+K_{33}^{(3)}& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 1, jj = 2
Total_K_Array(3, 4) = Total_K_Array(3, 4) + Element_K_Array_3(1, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}+K_{33}^{(3)}& K_{34}^{(3)}\\ 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 2, jj = 1
Total_K_Array(4, 3) = Total_K_Array(4, 3) + Element_K_Array_3(2, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}+K_{33}^{(3)}& K_{34}^{(3)}\\ 0& 0& K_{43}^{(3)}& 0\end{array}\right]\end{array}\end{array}$$

ii = 2, jj = 2
Total_K_Array(4, 4) = Total_K_Array(4, 4) + Element_K_Array_3(2, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}+K_{33}^{(3)}& K_{34}^{(3)}\\ 0& 0& K_{43}^{(3)}& K_{44}^{(3)}\end{array}\right]\end{array}\end{array}$$

最终得到整体刚度矩阵如上。组合方式如下图所示

图2 一维多连杆结构整体刚度矩阵组装方式

1.3 单元应力

Element_Stress = Bar1D2Node_Stress(Element_Dis_Vec, Geo_L)

• ToDo：计算单元的应力
• Input：单元位移列阵Element_Dis_Vec、长度Geo_L
• Output：单元应力标量Element_Stress

function Element_Stress = Bar1D2Node_Stress(Element_Dis_Vec, Geo_L)
    Element_Stress = [-1/Geo_L, 1/Geo_L] * Element_Dis_Vec
end

求取第1根杆件上的单元应力为

Element_Stress = Bar1D2Node_Stress(Element_Dis_Vec, Geo_L)

$$\begin{array}{c}\begin{array}{rl}{\sigma }^{(1)}& =\left[\begin{array}{cc}-1/l& 1/l\end{array}\right]\ast {\mathbf{q}}^{(1)}\\ & =\left[\begin{array}{cc}-1/l& 1/l\end{array}\right]\ast \left[\begin{array}{c}{u}_1^{(1)}\\ u_2^{(1)}\end{array}\right]\\ & ={\sigma }^{(1)}\end{array}\end{array}$$

1.4 单元节点力矢量

Element_Force_Array = Bar1D2Node_Force(Element_K_Array, Element_Dis_Vec)

• ToDo：计算单元的节点力矢量
• Input：单元刚度矩阵Element_K_Array、单元位移列阵Element_Dis_Vec
• Output：单元节点力矢量Element_Force_Array

function Element_Force_Array = Bar1D2Node_Force(Element_K_Array, Element_Dis_Vec)
    Element_Force_Array = Element_K_Array * Element_Dis_Vec
end

求取第1根杆件上的节点力为

Element_Force_Array = Bar1D2Node_Force(Element_K_Array_1, Element_Dis_Vec)

$$\begin{array}{c}\begin{array}{rl}{\mathbf{P}}^{(1)}& ={\mathbf{K}}^{(1)}\ast {\mathbf{q}}^{(1)}\\ & =\left[\begin{array}{cc}{K}_{11}^{(1)}& K_{12}^{(1)}\\ K_{21}^{(1)}& K_{22}^{(1)}\end{array}\right]\ast \left[\begin{array}{c}{u}_1^{(1)}\\ u_2^{(1)}\end{array}\right]\\ & =\left[\begin{array}{c}{P}_1^{(1)}\\ P_2^{(1)}\end{array}\right]\end{array}\end{array}$$

求得分别为第1根杆件上2个节点的节点力。

2 平面(2D)杆单元的有限元分析程序(Bar2D2Node)

2.1 单元刚度矩阵

Element_K_Array = Bar2D2Node_Stiffness(Mat_E, Geo_A, Node1_X, Node1_Y, Node2_X, Node2_Y)

• ToDo：计算单元的刚度矩阵
• Input：弹性模量Mat_E、横截面积Geo_A、节点1坐标(Node1_X, Node1_Y)、节点2坐标(Node2_X, Node2_Y)
• Output：单元刚度矩阵Element_K

function Element_K_Array = Bar2D2Node_Stiffness(Mat_E, Geo_A, Node1_X, Node1_Y, Node2_X, Node2_Y)
    % 杆件长度
    Geo_L = sqrt((Node2_X - Node1_X)^2 + (Node2_Y - Node1_Y)^2);
    % 杆件方向向量与整体坐标系X轴的夹角
    Alpha_Rad = atan2(Node2_Y - Node1_Y, Node2_X - Node1_X);
    Cos_Alpha = cos(Alpha_Rad);
    Sin_Alpha = sin(Alpha_Rad);
    Element_K_Array = Mat_E * Geo_A / Geo_L * ...
    [Cos_Alpha^2, Cos_Alpha * Sin_Alpha, -Cos_Alpha^2, -Cos_Alpha * Sin_Alpha; ...
    Cos_Alpha * Sin_Alpha, Sin_Alpha^2, -Cos_Alpha * Sin_Alpha, -Sin_Alpha^2; ...
    -Cos_Alpha^2, -Cos_Alpha * Sin_Alpha, Cos_Alpha^2, Cos_Alpha * Sin_Alpha; ...
    -Cos_Alpha * Sin_Alpha, -Sin_Alpha^2, Cos_Alpha * Sin_Alpha, Sin_Alpha^2];
end

图3 平面三连杆结构

上图为平面三连杆结构，通过程序得到每个杆件的单元刚度方程分别为

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{(1)}=\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}& K_{34}^{(1)}\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}& K_{44}^{(1)}\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{(2)}=\left[\begin{array}{cccc}{K}_{33}^{(2)}& K_{34}^{(2)}& K_{35}^{(2)}& K_{36}^{(2)}\\ K_{43}^{(2)}& K_{44}^{(2)}& K_{45}^{(2)}& K_{46}^{(2)}\\ K_{53}^{(2)}& K_{54}^{(2)}& K_{55}^{(2)}& K_{56}^{(2)}\\ K_{63}^{(2)}& K_{64}^{(2)}& K_{65}^{(2)}& K_{66}^{(2)}\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{(3)}=\left[\begin{array}{cccc}{K}_{55}^{(3)}& K_{56}^{(3)}& K_{57}^{(3)}& K_{58}^{(3)}\\ K_{65}^{(3)}& K_{66}^{(3)}& K_{67}^{(3)}& K_{68}^{(3)}\\ K_{75}^{(3)}& K_{76}^{(3)}& K_{77}^{(3)}& K_{78}^{(3)}\\ K_{85}^{(3)}& K_{86}^{(3)}& K_{87}^{(3)}& K_{88}^{(3)}\end{array}\right]\end{array}\end{array}$$

2.2 整体刚度矩阵

Total_K_Array = Bar2D2Node_Assembly(Total_K_Array, Element_K_Array, ElementNode_1, ElementNode_2)

• ToDo：计算整体刚度矩阵
• Input：单元刚度矩阵Element_K_Array、单元两端节点编号ElementNode_1、ElementNode_2
• Output：单元刚度矩阵Total_K_Array

function Total_K_Array = Bar2D2Node_Assembly(Total_K_Array, Element_K_Array, ElementNode_1, ElementNode_2)
    Element_DOF_Array(1) = 2 * ElementNode_1 - 1;
    Element_DOF_Array(2) = 2 * ElementNode_1;
    Element_DOF_Array(3) = 2 * ElementNode_2 - 1;
    Element_DOF_Array(4) = 2 * ElementNode_2;
    for ii = 1 : 4
        for jj = 1 : 4
            Total_K_Array(Element_DOF_Array(ii), Element_DOF_Array(jj)) = Total_K_Array(Element_DOF_Array(ii), Element_DOF_Array(jj)) + Element_K_Array(ii, jj);
        end
    end
    Total_K_Array = Total_K_Array;
end

2.2.1 组装第1根杆件

通过程序，初始时Total_K_Array为$8\times 8$的空矩阵，当组装第一根杆件时，有

Total_K_Array = Bar2D2Node_Assembly(Total_K_Array, Element_K_Array_1, 1, 2);

此时，

Element_DOF_Array(1) = 1;
Element_DOF_Array(2) = 2;
Element_DOF_Array(3) = 3;
Element_DOF_Array(4) = 4;

循环有，

ii = 1, jj = 1
Total_K_Array(1, 1) = Total_K_Array(1, 1) + Element_K_Array_1(1, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 1, jj = 2
Total_K_Array(1, 2) = Total_K_Array(1, 2) + Element_K_Array_1(1, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

组装完第1根杆后得到

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}& K_{34}^{(1)}& 0& 0& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}& K_{44}^{(1)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

2.2.2 组装第2根杆件

当组装第2根杆件时，有

Total_K_Array = Bar2D2Node_Assembly(Total_K_Array, Element_K_Array_2, 2, 3);

此时，

Element_DOF_Array(1) = 3;
Element_DOF_Array(2) = 4;
Element_DOF_Array(3) = 5;
Element_DOF_Array(4) = 6;

循环有，

ii = 1, jj = 1
Total_K_Array(3, 3) = Total_K_Array(3, 3) + Element_K_Array_2(1, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}+K_{33}^{(2)}& K_{34}^{(1)}& 0& 0& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}& K_{44}^{(1)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 1, jj = 2
Total_K_Array(3, 4) = Total_K_Array(3, 4) + Element_K_Array_2(1, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}+K_{33}^{(2)}& K_{34}^{(1)}+K_{34}^{(2)}& 0& 0& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}& K_{44}^{(1)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

组装完第2根杆后得到

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}+K_{33}^{(2)}& K_{34}^{(1)}+K_{34}^{(2)}& K_{35}^{(2)}& K_{36}^{(2)}& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}+K_{43}^{(2)}& K_{44}^{(1)}+K_{44}^{(2)}& K_{45}^{(2)}& K_{46}^{(2)}& 0& 0\\ 0& 0& K_{53}^{(2)}& K_{54}^{(2)}& K_{55}^{(2)}& K_{56}^{(2)}& 0& 0\\ 0& 0& K_{63}^{(2)}& K_{64}^{(2)}& K_{65}^{(2)}& K_{66}^{(2)}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

2.2.3 组装第3根杆件

当组装第3根杆件时，有

Total_K_Array = Bar2D2Node_Assembly(Total_K_Array, Element_K_Array_3, 3, 4);

此时，

Element_DOF_Array(1) = 5;
Element_DOF_Array(2) = 6;
Element_DOF_Array(3) = 7;
Element_DOF_Array(4) = 8;

循环有，

ii = 1, jj = 1
Total_K_Array(5, 5) = Total_K_Array(5, 5) + Element_K_Array_2(1, 1);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}+K_{33}^{(2)}& K_{34}^{(1)}+K_{34}^{(2)}& K_{35}^{(2)}& K_{36}^{(2)}& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}+K_{43}^{(2)}& K_{44}^{(1)}+K_{44}^{(2)}& K_{45}^{(2)}& K_{46}^{(2)}& 0& 0\\ 0& 0& K_{53}^{(2)}& K_{54}^{(2)}& K_{55}^{(2)}+K_{55}^{(3)}& K_{56}^{(2)}& 0& 0\\ 0& 0& K_{63}^{(2)}& K_{64}^{(2)}& K_{65}^{(2)}& K_{66}^{(2)}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

ii = 1, jj = 2
Total_K_Array(5, 6) = Total_K_Array(5, 6) + Element_K_Array_2(1, 2);

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}+K_{33}^{(2)}& K_{34}^{(1)}+K_{34}^{(2)}& K_{35}^{(2)}& K_{36}^{(2)}& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}+K_{43}^{(2)}& K_{44}^{(1)}+K_{44}^{(2)}& K_{45}^{(2)}& K_{46}^{(2)}& 0& 0\\ 0& 0& K_{53}^{(2)}& K_{54}^{(2)}& K_{55}^{(2)}+K_{55}^{(3)}& K_{56}^{(2)}+K_{56}^{(3)}& 0& 0\\ 0& 0& K_{63}^{(2)}& K_{64}^{(2)}& K_{65}^{(2)}& K_{66}^{(2)}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$

组装完第3根杆后得到

$$\begin{array}{c}\begin{array}{r}{\mathbf{K}}^{All}=\left[\begin{array}{cccccccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}& 0& 0& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}& 0& 0& 0& 0\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}+K_{33}^{(2)}& K_{34}^{(1)}+K_{34}^{(2)}& K_{35}^{(2)}& K_{36}^{(2)}& 0& 0\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}+K_{43}^{(2)}& K_{44}^{(1)}+K_{44}^{(2)}& K_{45}^{(2)}& K_{46}^{(2)}& 0& 0\\ 0& 0& K_{53}^{(2)}& K_{54}^{(2)}& K_{55}^{(2)}+K_{55}^{(3)}& K_{56}^{(2)}+K_{56}^{(3)}& K_{57}^{(3)}& K_{58}^{(3)}\\ 0& 0& K_{63}^{(2)}& K_{64}^{(2)}& K_{65}^{(2)}+K_{65}^{(3)}& K_{66}^{(2)}+K_{66}^{(3)}& K_{67}^{(3)}& K_{68}^{(3)}\\ 0& 0& 0& 0& K_{75}^{(3)}& K_{76}^{(3)}& K_{77}^{(3)}& K_{78}^{(3)}\\ 0& 0& 0& 0& K_{85}^{(3)}& K_{86}^{(3)}& K_{87}^{(3)}& K_{88}^{(3)}\end{array}\right]\end{array}\end{array}$$

最终得到整体刚度矩阵如上。组合方式如下图所示

图4 平面多连杆结构整体刚度矩阵组装方式

2.3 单元应力

Element_Stress = Bar2D2Node_Stress(Mat_E, Node1_X, Node1_Y, Node2_X, Node2_Y, Element_Dis_Vec)

• ToDo：计算单元的应力
• Input：结构弹性模量Mat_E、第1个节点坐标Node1_X, Node1_Y、第2个节点坐标Node2_X, Node2_Y、单元位移列阵Element_Dis_Vec
• Output：单元应力标量Element_Stress

function Element_Stress = Bar2D2Node_Stress(Mat_E, Node1_X, Node1_Y, Node2_X, Node2_Y, Element_Dis_Vec)
    % 杆件长度
    Geo_L = sqrt((Node2_X - Node1_X)^2 + (Node2_Y - Node1_Y)^2);
    % 杆件方向向量与整体坐标系X轴的夹角
    Alpha_Rad = atan2(Node2_Y - Node1_Y, Node2_X - Node1_X);
    Cos_Alpha = cos(Alpha_Rad);
    Sin_Alpha = sin(Alpha_Rad);
    Trans_Array = [Cos_Alpha, Sin_Alpha, 0, 0; 0, 0, Cos_Alpha, Sin_Alpha];
    Element_Stress = [-1/Geo_L, 1/Geo_L] * Trans_Array * Element_Dis_Vec;
end

求取第1根杆件上的单元应力为

Element_Stress = Bar2D2Node_Stress(Mat_E, Node1_X, Node1_Y, Node2_X, Node2_Y, Element_Dis_Vec)

$$\begin{array}{c}\begin{array}{rl}{\sigma }^{(1)}& =\left[\begin{array}{cc}-1/l& 1/l\end{array}\right]\ast {\mathbf{T}}^{(1)}\ast {\overline{\mathbf{q}}}^{(1)}\\ & =\left[\begin{array}{cc}-1/l& 1/l\end{array}\right]\ast \left[\begin{array}{cccc}\cos \alpha & \sin \alpha & 0& 0\\ 0& 0& \cos \alpha & \sin \alpha \end{array}\right]\ast \left[\begin{array}{c}{u}_1^{(1)}\\ v_1^{(1)}\\ u_2^{(1)}\\ v_2^{(1)}\end{array}\right]\\ & ={\sigma }^{(1)}\end{array}\end{array}$$

2.4 单元节点力矢量

Element_Force_Array = Bar2D2Node_Force(Element_K_Array, Element_Dis_Vec)

• ToDo：计算单元的节点力矢量
• Input：单元刚度矩阵Element_K_Array、单元位移列阵Element_Dis_Vec
• Output：单元应力节点力矢量Element_Force_Array

function Element_Force_Array = Bar2D2Node_Force(Element_K_Array, Element_Dis_Vec)
    Element_Force_Array = Element_K_Array * Element_Dis_Vec
end

求取第1根杆件上的节点力为

Element_Force_Array = Bar2D2Node_Force(Element_K_Array_1, Element_Dis_Vec)

$$\begin{array}{c}\begin{array}{rl}{\mathbf{P}}^{(1)}& ={\mathbf{K}}^{(1)}\ast {\mathbf{q}}^{(1)}\\ & =\left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& K_{13}^{(1)}& K_{14}^{(1)}\\ K_{21}^{(1)}& K_{22}^{(1)}& K_{23}^{(1)}& K_{24}^{(1)}\\ K_{31}^{(1)}& K_{32}^{(1)}& K_{33}^{(1)}& K_{34}^{(1)}\\ K_{41}^{(1)}& K_{42}^{(1)}& K_{43}^{(1)}& K_{44}^{(1)}\end{array}\right]\ast \left[\begin{array}{c}{u}_1^{(1)}\\ v_1^{(1)}\\ u_2^{(1)}\\ v_2^{(1)}\end{array}\right]\\ & =\left[\begin{array}{c}{P}_{1x}^{(1)}\\ P_{1y}^{(1)}\\ P_{2x}^{(1)}\\ P_{2y}^{(1)}\end{array}\right]\end{array}\end{array}$$

求得分别为第1根杆件上2个节点的节点力。

参考

1. 曾攀. 有限元分析基础教程[M]. 北京: 清华大学出版社, 2008.