有限元进阶编程——预处理法组装整体刚度矩阵

前言

我们的有限元基础编程课告一段落，接下来将会为大家带来有限元进阶编程，采用的编程语言依然是Matlab，精彩多多，感兴趣的同学可以持续关注哦~

本次推文给大家分享的是整体刚度矩阵组装另一种方法——预处理法，该方法采用定位向量组装形成整体刚度矩阵，规模较小，考虑了支座约束的 0 位移自由度，，求解未知位移时，可直接对刚度矩阵进行线性方程组求解，避免了边界条件特殊处理（划 0 置 1 、乘大数法）。

预处理法

确定单元刚度矩阵各元素在总刚矩阵的位置，就可以由单元刚度矩阵直接集成总刚矩阵。

定位向量的定义：单元杆端位移分量（局部位移码）所对应的结构节点位移分量（整体 位移码）序号组成的向量，称为单元的定位向量。概念大概看一下，具体看下面操作

例题：

如上图所示，结构各单元的定位向量为（被约束的位移编码一般记为0）

单元 1（0,0,0,1,2,3）；单元 2（1,2,3,4,5,6）；单元 3（4,5,6,0,0,0）

具体操作：

1. 求出单元刚度矩阵；

2. 将单元的定位向量分别写在单元刚度矩阵的上方和右侧，使得单元刚度矩阵的分量与单元定位向量的分量一一对应；

   

   
   

   

3. 单元定位向量为 0 所对应的分量，在单元刚度矩阵相应的行和列可以删掉；

4. 对于单元定位向量中不为 0 的分量，对应的就是单元刚度矩阵在整体刚度矩阵中的行码和列码。 

   
   

   

相关代码： 

1. 先将边界条件矩阵初始化为单位矩阵，然后将支座约束节点位移设置为 0 ，最后将自由节点位移信息按照 1，2，3...顺序来编排。

   
   

2.  变量：nf(边界条件矩阵);nnd(节点数目);nodof(节点自由度数目);node_1(单元第一个节点)

nf = ones(nnd, nodof); % Initialise the matrix nf to 1
nf(1,2) =0 ;               % Prescribed nodal freedom of node 1
nf(6,1)= 0 ; nf(6,2)= 0 ;  % Prescribed nodal freedom of node 6
% Counting of the free degrees of freedom
n=0;
for i=1:nnd
    for j=1:nodof
        if nf(i,j) ~= 0
        n=n+1;
        nf(i,j)=n;
        end
    end
end
g=[nf(node_1,1); nf(node_1,2); nf(node_2,1);nf(node_2,2)]

以上代码取自《Matlab与Abaqus有限元分析理论与应用》（汉译版），代码风格的源头是《Programming the Finite Element Method》，特别说明的是：该书是有限元编程的“里程碑”著作，内容非常全面，代码量巨大，免费开源，但是对于我们新手很不友好，代码语言采用Fortran编写，阅读起来十分费劲，后续不断有人在该著作的基础上将代码语言修改为 Matlab 语言，与 Fortran 相比起来，易读性显著提升，可视化更强。好啦好啦，题外话不多说啦，继续我们的正题！

总结：利用定位向量组装整体刚度矩阵，考虑了支座约束的节点 0 自由度，在组装之前就把含有 0 位移的分量“划掉”，组装后的整体刚度矩阵规模较小，对应的位移列阵皆为未知位移，等效节点载荷列阵也皆为已知量，在求解位移时，直接对线性方程组求解即可，避免了边界条件特殊处理（划 0 置 1 、乘大数法）。

从以上描述可知，该方法又名为“预处理法”的缘故。更加高阶的技巧：整体刚度变带宽存储、线性方程组解法等将在以后为大家呈现~

-End-