ANSYS模型刚度、质量矩阵快速提取小软件—km_from_Ansys

1. 背景

从事结构振动控制、车桥耦合振动、结构健康监测传感器优化布置、结构动力性能分析等等一系列研究的同仁们应该都面临过一个同样的问题—“怎么把结构的刚度和质量矩阵建立出来？”。这对于那些数值分析高手和专家可能不是什么问题；但是对于科研刚入门的新手来说，这个难度还是相当大的。如果都靠自己写程序来建立有限元模型，则对理论基础、编程水平都有很高的要求，甚至程序做出来也未必能保证其正确性，是一个很让人头疼的问题。

对于一些简单的被动控制装置或简单的动力学分析，当然也可以在有限元分析软件中构造出装置组成直接分析（刚度+阻尼类型），但是对于稍复杂一些的控制装置和耦合分析等问题，会受到平台功能上的限值，尤其是对于主动和半主动等涉及控制算法的研究来说，基本很难在有限元软件平台上实现分析。再加上如果需要对装置进行参数优化，需要进行多次重复计算，难度就更大。

Ansys、ABAQUS等软件平台给我们提供了比较稳定有效的有限元模型建立平台，通过借助商业软件来建立模型，再将其中的刚度、质量矩阵导出，是非常可取的一种方法。如果能够提取出模型的矩阵，明晰计算原理，就能够很容易的通过自己的程序设计对计算过程进行补充、调整，来达到自己定制的计算分析目的。其实，不仅对于振动控制，比如结构静动力分析、车桥耦合分析、结构传感器优化配置方案设计等，都有应用需求。因此，一个能够便捷的提取结构矩阵的方法就显得至关重要。

技术邻平台已经有大佬提供了ABAQUS软件刚度和质量矩阵的导出方法。这里补充一下在ANSYS中导出质量和刚度矩阵的方法和小软件。

2. 提取方法

Ansys中其实提供了几种单元刚度矩阵和整体刚度矩阵的提取接口，先做一下简单介绍：

（1）单元刚度矩阵提取，集成整体刚度矩阵

首先选择要提取的单元，通过/debug命令输出选择集中的单元刚度矩阵，再根据结构的单元连接信息，组合成整体刚度矩阵，命令如下：

/OUTPUT,cp,out,, ! 将输出信息送到cp.out文件 

/debug,-1,,,1 ! 指定输出单元矩阵 

/SOLU 

SOLVE 

finish 

/OUTPUT, TERM ! 将输出信息送到output windows中 

感受：这个方法很直接，但稍显麻烦，需要提取到结构的单元组成信息、节点信息以及单元刚度矩阵等，对于稍大一些结构来说，可能需要1个小时甚至更长时间才能完成。

（2）超单元方法

/solu   

antype,7 !substructuring分析类型   

seopt,matname,1 !设置文件名称和刚度矩阵类型(刚度，质量，阻尼等)   

nsel,all !选择所有节点   

m,all,all !定义所有节点自由度为主自由度   

solve !求解   

selist,matname,3 !列出整体刚度矩阵  

这种方法是可以提取到结构整体刚度矩阵的，但是问题在于需要指定主自由度，对于一般结构而言，这个方法还是适用的，问题是如果结构中存在耦合关系，结构自由度存在依赖关系，该方法就不太适用了，可能还有其他解决办法，但是我后来还是放弃这种方法。

（3）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中的矩阵文件是以稀疏矩阵的形式存储的，得到的十进制文件，只是给出了矩阵中元素对应的位置和值，我们需要对这一串数字进行整理才能得到原始的矩阵形式。这就需要明确十进制文件中数字的排列规律。

文件基本格式是前面有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格式则仅记录下三角元素。

        根据Harwell-Boeing文件格式，可读取矩阵的任意行列元素的数值，也可编程还原为满矩阵存储，以便它用，很显然这种提取方式比较方便。如当生成.FULL文件后，可采用命令/AUX2$FILE,mywork,full$HBMAT,mystiff,txt,ASCII,STIFF,YES$FINISH将二进制mywork.full文件输出为ASCII码文件mystiff.txt，并输出右边项。

按照这个规律，对十进制文件中的数字进行整理，就可以得到想要的矩阵了。说到这里其实已经结束了，但这只是方法，从方法到应用还是需要一定的摸索和尝试，为了便于应用，下面介绍一个可以便捷提取的小工具。

3. 快速提取小软件

该工具是基于上述第三种方法做的简单开发，目的是提高工作效率，把时间更多的用在需要做的理论研究上面，所以这个工具会尽可能的减少复杂操作，尽可能傻瓜式，尽量使User无需看上面的那些繁复介绍就可以实现矩阵的提取。工具基于Matlab软件开发，长相一般，是这个样......,

但是功能是：五个按键解决提取问题。

4. 应用实例

(1). 首先在ansys中建立一个有限元模型，这里以一个简单3层框架为例：

fini

/clear

/filename,'Frame'          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    文件名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

运行完成后，在根目录下生成了四个文件，这便是转换完成的十进制文件，接着需要做的就是将十进制文件转换为矩阵形式，就是通过工具来实现。

(3). 通过工具实现提取

A. 打开运行工具，在km_from_Ansys程序目录下，在Matlab命令行输入km_from_Ansys，是初始界面

B.点击“刚度文件路径”后的“获取路径按钮”，找到刚刚生成的十进制文件的存储路径，选择刚度十进制文件Stiff.txt：

同样的方法选择好质量矩阵的文件路径，和Mapping文件路径（Mapping文件有两个，选择任意一个都可以），以及想要把提取完成后的矩阵文件存储的路径，全部选择完成后，点击“生成文件”。

C. 接着是大约几秒的工具处理过程，根据模型大小不同，所需时间不等。处理完成后，根据提示输入生成的矩阵文件的存储文件名，比如“KM”。

，点击“OK”，完成操作，提示信息为文件的存储路径，文件格式为Matlab数据文件格式。

D. 验证

Matlab打开生成的矩阵文件，一共有三个矩阵，分别是kk,mm,address,其中kk文件和mm文件分别是结构刚度和质量矩阵文件，以稀疏矩阵形式存储，address文件是矩阵的自由度说明，其中第二列是对应节点号，第三列是对应自由度编码，1,2,3,4,5,6,分别代表X,Y,Z,ROTX,ROTY和ROTZ。便于对矩阵进行操作。

通过求取模态特征值来验证提取矩阵的有效性：

>>load('KM.mat')

>> [v,d]=eigs(kk,mm,10,'sm');

>> d=sqrt(d);

>> d=d/2/pi;

>> d=sum(d);

>> d=d';

上图为Ansys中模态分析的结果，下图为提取矩阵的分析结果，结果一致，表明提取矩阵有效。

6. 说明

提取工具适用于任意单元，任意结构形式，但仅限于弹性模型。