有限元进阶编程——一维变带宽存储整体刚度矩阵

1年前浏览721

编辑：易木木响叮当

一维变带宽存储整体刚度矩阵

必要性探讨

在有限元分析过程中，整体刚度的组装是必不可少的步骤，所谓的“组装”就是将单元刚度矩阵各个分量按照一定的顺序放在整体刚度矩阵中，由于仅有在同一单元中相关联的节点才会在整体矩阵相应的行及列中出现刚度系数，则在整体矩阵中，必然会有大量的零数据。若以全矩阵的形式存储刚度矩阵 K ，随着求解规模的增大，单元数目势必会增多，将会导致内存占用量急剧增大。以四节点四边形单元（每个节点两个自由度）为例进行说明。

假设模型有 n 个单元，则节点数为：

，整体刚度矩阵阶数为：

，则保存在内存中的总数居个数为：

。

数据若以双精度浮点数的形式存储，则每个数据占用8字节(Byte)，那么总的内存消耗量 Memory 为

根据上式，可绘制以下曲线图：

图1 整体刚度矩阵存储量

由上图可知，当单元数目增加到 2 万个时，整体刚度矩阵的内存占用量已经超过 12 GB，因此有限元方法应用于大规模求解模型时，需将整体刚度的存储方式更改为稀疏矩阵存储、或者按照刚度矩阵带宽的特点进行存储。但是目前的Matlab中，尚不能指定稀疏矩阵是对称的，会浪费不少内存。因此本次推文将介绍存储量最小的方式：一维变带宽紧缩存储。

一维变带宽

带宽

我们先来看一下带宽的概念，下图为一平面问题的三角形单元离散，第

个单元的节点位移列阵为：

该单元以全矩阵的形式装配时在整体刚度矩阵中对应的位置：

图2 单元刚度矩阵带宽表示方式

从以上图可以大致了解单元刚度矩阵带宽的概念。

一维变带宽存储

一维变带宽存储的存储方式：

图3 一维变带宽存储方法示意图

注：

1. 只存储对称矩阵的下三角区（或上三角区）的元素；
2. 整体刚度矩阵下三角区的每一行（或每一列）第一个非零元素至对角线元素，按序排列在一维数组 kv 当中；
3. 将主对角线的元素在数组 kv 中的位置存放在一维数组 kdiag 中，K矩阵中第行第列元素（）在数组 kv 中的地址 :
以上就是一维变带宽存储的理论概念，接下来我们将借助《Programming the finite element method》中的平面桁架案例来比较全矩阵存储、预处理刚度矩阵存储、一维变带宽存储三种存储刚度矩阵的方法。

模型

图4 《Programming the finite element method》 P41

整体刚度矩阵大小对比

直接刚度法	12*12 double
预处理法	9*9 double
一维变带宽+预处理法	1*39 double

图5 直接刚度法存储整体刚度矩阵

图6 预处理法法存储整体刚度矩阵

图7 一维变带宽+预处理法存储整体刚度矩阵

代码

代码取自《有限单元法基础及Matlab编程第二版》

% --------一维变带宽存储整体刚度矩阵------------
function kv = fsparv(kv,km,g,kdiag)
% FSPARV 函数功能：
%   集装总体刚度矩阵
% 输入:
%   kdiag 主对角元素地址向量
%   km 单元刚度矩阵
%   g 单元每个结点的自由度编号，也即定位向量
%   kdiag 主对角元素地址向量
% 输出:
%   kv 总体刚度矩阵
% -----------------------------------------------------
    idof = length(g);          % 定位向量尺寸
    for i = 1 : idof
        k = g(i);
        if k ~= 0              % 定位向量元素非零
            for j = 1 : idof
                if g(j) ~= 0   % 定位向量元素非零
                    iw = k - g(j);  % 自由度号差值
                    if iw >= 0
                        ival = kdiag(k) - iw;
                        kv(ival) = kv(ival) + km(i,j);
                    end
                end
            end
        end
    end
end

谢谢你看完木木同学的分享，今日份阅读花费的流量+1M哈哈哈哈哈哈。

-End-

易木木响叮当

想陪你一起度过短暂且漫长的科研生活

来源：易木木响叮当

理论

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2023-06-01