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. 說明

提取工具適用于任意單元，任意結構形式，但僅限于彈性模型。