提取ANSYS中的刚度矩阵和质量矩阵到MATLAB中进行二次开发

最近在考虑自己编写的程序和商用软件的验证问题，有限元结构分析中最关键的一环就是刚度矩阵的获得，如果涉及到模态分析，还有质量矩阵。考虑到商业软件的成熟性，可以用ANSYS生成的刚度矩阵做参照来看自己编写的程序是否正确，因此如何提取ANSYS中结构的刚度矩阵，并进行随后的验证或者二次开发是一个问题。

https://www.jishulink.com/content/post/1796144

受上述帖子的启发，使用MATLAB提取ANSYS中的刚度和质量矩阵，并进行模态分析验证提取的矩阵的正确性。

首先，在ANSYS中使用HBMAT命令方法提取整体矩阵。

命令：HBMAT,fname,ext,--,form,matrx,rhs

其中：

Fname---输出矩阵的路径和文件名，缺省为当前工作路径和当前工作文件名。

ext---输出矩阵文件的扩展名，缺省为.matrix。

form---定义输出矩阵文件的格式，其值可取：

=ASCII：ASCII码格式；

=BIN：二进制格式。

matrix---定义输出矩阵的类型，其值可取：

=STIFF：输出刚度矩阵。可用于写入了.FULL文件的任何类型的分析。

=MASS：输出质量矩阵。可用于特征值屈曲、子结构分析、模态分析。

=DAMP：输出阻尼矩阵。仅用于有阻尼的模态分析。

rhs---右边项输出控制（右边项指用矩阵所表示方程的等号右端矢量，这里可为节点荷载向量），如rhs=YES则输出，如rhs=NO则不输出。

模态分析时，因仅LANB和QR法可生成完整的质量矩阵，因此也仅采用这两种方法时才可使用HBMAT命令得到质量矩阵文件。

这个是Ansys帮助文件中对HBMAT命令的解释，在Ansys分析中，会在根目录下形成.full的二进制文件，里面存储的就是已经组合好的刚度和质量矩阵，只是是以二进制文件的形式存在，通过HBMAT命令可以将二进制文件转换为可以使用的十进制文件，十进制文件是以txt的形式存放在根目录下。

ANSYS中根据Harwell-Boeing文件格式存储刚度矩阵，如图所示：

文件基本格式是前面有4或5行描述数据，其后为单列矩阵元素值，说明如下：

第1行：格式（A72），为文件头的字符型解释，如刚度矩阵或质量矩阵等标题。

第2行：格式（5I14），分别表示该文件的总行数（不包括文件头）、矩阵列指针的总行数、矩阵行索引的总行数、矩阵元素数值的总行数、右边项总行数。

第3行：格式（A3,11X,4I14），分别为矩阵类型、矩阵行数、矩阵列数、矩阵行索引数（对组装后的矩阵，该值等于矩阵行索引数）、单元元素数（对组装后的矩阵此值为0）。

第4行：格式（2A16,2A20），分别表示列指针格式、行索引格式、系数矩阵数值格式、右边项数值格式。

第5行：格式（A3,11X,2I14），A3各列分别表示右边项格式、应用高斯起始矢量、应用eXact求解矢量；两个整数分别表示右边项列数、行索引数。三个字符中的第1个字符可取：F---全部存贮（如节点荷载向量的全部元素）、M---与系数矩阵相同方法。

第6行后：矩阵元素值（单列）。

矩阵类型用3个字符表示，第1个字符可取：R---实数矩阵、C---复数矩阵、P---仅矩阵结构（无元素数值）；第2个字符可取：S---对称矩阵、U---不对称矩阵、H---Hermitian矩阵、Z---病态对称矩阵；R---带状矩阵；第3个字符可取：A---组装的矩阵、E---单元矩阵（未组装）。对称矩阵只存储下三角元素，如结构刚度矩阵为对称矩阵，Harwell-Boeing格式则仅记录下三角元素。

通过上述解释我们可以知道该刚度矩阵有1027个列指针，有7978个行索引，有1026个右边项，矩阵的维度为1026×1026。

要将上述文件读取到MATLAB中并转化为常规的稀疏矩阵，需要了解Harwell-Boeing压缩储存格式及重构方法，这里面采用的是CSC方式，以列指针、行索引的方式去重构稀疏矩阵。

以三层框架结构为参考例子，其模态分析命令流如下：

fini

/clear

/filename,'Frame'    

/prep7

et,1,beam188

et,2,shell63

mp,ex,1,3.0e10

mp,nuxy,1,0.2

mp,dens,1,2500

r,1,0.15,0.15,0.15,0.15

sectype,1,beam,rect

secdata,0.8,0.8

k,1,

k,2,,,5

k,3,10,,5

k,4,10

k,5,,3.5

k,6,,3.5,5

k,7,10,3.5,5

k,8,10,3.5

k,9,5,3.5,,

k,10,5,3.5,5

k,11,,7

k,12,,7,5

k,13,10,7,5

k,14,10,7

k,15,5,7,,

k,16,5,7,5

k,17,,10.5

k,18,,10.5,5

k,19,10,10.5,5

k,20,10,10.5

k,21,5,10.5,,

k,22,5,10.5,5

l,1,5

l,2,6

l,3,7

l,4,8

l,5,11

l,6,12

l,7,13

l,8,14

l,11,17

l,12,18

l,13,19

l,14,20

a,5,6,10,9

a,7,8,9,10

a,11,12,16,15

a,13,14,15,16

a,17,18,22,21

a,19,20,21,22

lsel,all

latt,1,,1,,,,1

lesize,all,,,4

lmesh,all

asel,all

aatt,1,1,2

esize,,4

amesh,all

allsel

nummrg,all

numcmp,node

nsel,s,loc,y,0

d,all,all

fini

!模态分析

/solu

allsel

antype,modal

modopt,lanb,10

solve

finish

!(2). 模态分析完成后，在ansys中运行HBMAT命令得到十进制文件：

! 获取矩阵十进制文件

/aux2

file,Frame,full !与文件名一致

hbmat,Stiff,txt,,ascii,stiffness,yes,yes   !提取刚度矩阵 存储为Stiff.txt

hbmat,Mass,txt,,ascii,mass,yes,yes     !提取质量矩阵  存储为Mass.txt

Finish

提取其前10阶模态如下：

通过MATLAB对文件提取并处理后，可以使用spy命令查看稀疏矩阵的分布，刚度矩阵分布和各变量维度如下：

稀疏刚度矩阵的具体分布如下：

稀疏质量矩阵的具体分布如下：

计算矩阵的特征值，利用提取的刚度矩阵和质量矩阵做模态分析。

与ANSYS中的结果一致，因此验证了正确性。

在Stiff.mapping文件中，格式是这样的：

在MATLAB中使用相同的方法提取并整理，如下：