提取刚度矩阵.docx

============为什么折腾这个文档========
我有一个计算线性动力学方程组的瞬态、谐响应和静力学的python程序，现希望开发一个将ANSYS组集好的总体矩阵导入该PYTHON程序中的接口。
该问题可分解为：
[STEP1]  [ANSYS]->[包含矩阵信息的文件]
[STEP2]  [包含矩阵信息的文件]->[python通用数据对象]  
[STEP3]  [python通用数据对象]->[程序特定数据对象]->[进行计算]
因此检索了一些帖子，基本上完成了这项工作，本文是对[STEP1]和[STEP2]的整理，并且利用[STEP3]对结果进行了验证
============主要内容==================
1，了解从ANSYS中提取总体矩阵和载荷向量的方法；
2，了解提取出来的矩阵是怎样表示的；
3，说明在Python中，如何读取这样的矩阵；
4，构造一个简单的算例，说明整个【建模】-【提取】-【读取】过程及其正确性；
=========站内检索综述====================
检索词：提取 矩阵
得到21个结果，代表性的帖子有下面这9个：

编号[1]
标题：ansys中怎样提取质量，刚度，阻尼矩阵？
地址： ... fromuid-159019.html
要点：pengweicai给出了一段网上最常见的提取代码,该程序以fortran写成，可以利用.full文件以及一些列约定将ANSYS中的总体矩阵读入FORTRAN中。
编号[2]
标题：如何得知HBMAT命令提取的质量、刚度矩阵对应的自由度？
地址： ... fromuid-159019.html
要点：提出了使用HBMAT命令提取稀疏矩阵时常见的问题：我们如何知道提取出来的信息是怎么储存的呢？
编号[3]
标题：[分享]ANSYS中整体、单元刚度和质量矩阵的提取
地址： ... fromuid-159019.html
要点：在该帖子的7楼，其实已经给出了帖子[2]中问题的解答，即HBMAT中提取出来的矩阵是Harwell-Boeing格式的，并且给出了该格式的细节，可惜是英文的，没引起多少关注。
编号[4]
标题：帮我看看提取的刚度与质量矩阵
地址： ... fromuid-159019.html
要点：这个帖子所示的矩阵并非是使用HBMAT命令提出出来的，而应该是SELIST命令列举出来的未压缩的矩阵，后续楼层的回帖给了大家一个提示，即有可能提取出来的矩阵是引入了边界条件的（即删除了被约束的行和列的）。
编号[5]
标题：提取刚度矩阵的问题
地址： ... fromuid-159019.html
要点：本帖作者的工作是基于单元刚度矩阵的，因此ANSYS中提取的单元刚度矩阵是否处于总体坐标系就成为问题。该问题并非本文内容，但仍值得关注。
编号[6]
标题：提取刚度矩阵丢失节点的问题
地址： ... fromuid-159019.html
要点：帖子[5]作者的又一帖，在这里帖子[5]的问题得到了欧阳中华老师的回答。
编号[7]
标题：提取刚度矩阵的ANSYS操作过程
地址： ... fromuid-159019.html
要点：实际上这就是使用HBMAT从ANSYS中提取总体矩阵的全过程！只是还有一些细节待确定。
编号[8]
标题：提取整体刚度矩阵、质量矩阵及阻尼矩阵的简单方法
地址： ... fromuid-159019.html
要点：给出了利用“不减缩的”子结构方法来得到总体矩阵的方法（这也是网络上常见的代码之一）
编号[9]
标题：质量矩阵、刚度矩阵如何提取？
地址： ... fromuid-159019.html
要点：16443在5楼的回帖中给出了提取刚度矩阵的三种方法

=======站外检索略述========================

百度检索：提取 矩阵
比较好的帖子有：
编号[10]
来源:百度文库
标题:怎样从ansys中提取单元刚度矩阵与质量矩阵
地址:
要点:这应该就是16443在帖子[9]中回复的内容了，全面的总结了在帖子[3,4,5,9]中涉及的问题。
编号[11]
来源:中华钢结构
标题:ansys刚度矩阵Harwell-Boeing格式的具体含义讨论
地址:http://okok.org/forum/viewthread.php?tid=184007
要点:如题，后续楼层给出了一些将矩阵读入ANSYS的APDL（好不容易读出来，又读进去干嘛呢……）
编号[12]
来源:simwe
标题:关于ANSYS（质量、刚度、阻尼）矩阵Harwell-boeing格式数据的说明
地址:
要点:比[11]更透彻的HB格式说明!

=============================================================
=======1.从ANSYS中提取总体矩阵的方法=================================
=============================================================
1,用/DEBUG命令
2,子结构法
3,HBMAT
详见帖子[10]
PS.个人感觉HBMAT方法最靠谱，一是它的格式（Harwell-boeing）在很多场合都是通用的，二是BHMAT命令是文档化的、功能就是用来提取总体刚度矩阵的命令。因此，相比于子结构法的剑走偏锋，/DEBUG命令的繁复，HBMAT命令方法更“标准”一些，因此在后文只关注此方法。
=======2.BH格式的矩阵是如何表示的===================================
HBMAT命令并不是很复杂的命令，稍复杂的地方是采用该命令提取出来的矩阵是经过压缩的，称为Harwell-boeing格式，也叫Compressed Sparse Column格式。
其具体压缩和还原方式见帖子[3]（English）或[11][12](中文)
=======3.如何在Python中读入BH格式的矩阵===============================上文说过，Harwell-boeing格式，也叫Compressed Sparse Column格式，而Python.scipy中就有这样的稀疏矩阵：
1. class scipy.sparse.csc_matrix(arg1, shape=None, dtype=None, copy=False, dims=None, nzmax=None)

可以通过HB文件中直接读取的行标指针，行标和数据创建，例如：
1. >>> indptr = array([0,2,3,6])
2. >>> indices = array([0,2,2,0,1,2])
3. >>> data = array([1,2,3,4,5,6])
4. >>> csc_matrix( (data,indices,indptr), shape=(3,3) ).todense()
5. matrix([[1, 0, 4],
6.         [0, 0, 5],
7.         [2, 3, 6]])

对应的HB文件应为(*号部分表示并非本例关注的数据)：
1. Rainyboy Testing Matrix in BH format         
2.            ***            4           6           6            0
3. RRA                       **          **         **            0
4. (I14)           (I14)           (d25.15)            (d25.15)            
5. 0
6. 2
7. 3
8. 6
9. 0
10. 2
11. 2
12. 0
13. 1
14. 2
15. 1
16. 2
17. 3
18. 4
19. 5
20. 6
由文件头可知，indptr的长度为4，因此0,2,3,6就是indotr的内容
indices的长度为6，因此后续的0,2,2,0,1,2就是indices的内容
data的长度为6，因此后续的1,2,3,4,5,6就是data的内容
=======4.一个【建模】-【提取】-【读取】-【计算】的例子===============
【建模APDL】
1. FINISH
2. /CLEAR
3. /TITLE,CASE STUDY _BEAM _BEAM3 BY RAINYBOY
4. /PREP7
5. /ESHAPE,1                 !显示壳单元厚度
6. !**********************
7. !几何参数表
8. !**********************
9. *SET,L_HORI,0.1                                                !横梁的长度
10. *SET,TA,0.005                                                  !正方形截面的边长
11. *SET,MESHCOUNT,2                                        !每段的分网数
12. *SET,IZZ,TA*TA*TA*TA/12                !转动惯量
13. *SET,IYY,TA*TA*TA*TA/12                !转动惯量
14. !**********************
15. !材料参数表
16. !**********************
17. *SET,MEX,1.78E11                                        !弹性模量
18. *SET,MPRXY,0.3                                                !泊松比
19. *SET,MDENS,7850                                                !密度
20. !**********************
21. !相关设置
22. !**********************
23. MP,EX,1,MEX                                                    !设置材料弹性模量
24. MP,PRXY,1,MPRXY                                                !设置材料泊松比
25. MP,DENS,1,MDENS                                                !设置材料密度
26. BETAD,1E-5                                                                !BETA阻尼
27. ET,1,BEAM3                                                    !设置平面梁单元
28. R,1,TA*TA,IZZ,TA                                        !设置截面参数
29. !DMPRAT,0.10000                                                !阻尼比
30. !**********************
31. !几何->分网
32. !**********************
33. TYPE,1                                                         !指定分网类型
34. MAT,1                                                                !指定材料类型
35. REAL,1                                                         !指定实参数
36. K,1,0,0,0                                                      !建立三个关键点
37. K,2,L_HORI,0,0                                                                                
38. L,1,2                                                                !建立几何体
39. ALLSEL,ALL
40. LESIZE,ALL,,,MESHCOUNT          !设置线段分网数
41. LMESH,ALL                                                      !分网
42. !**********************
43. !几何约束
44. !**********************
45. ALLSEL,ALL
46. NSEL,S,LOC,X,0                                                !选择固定端节点
47. D,ALL,ALL                                                      !设置为约束所有自由度
48. ALLSEL,ALL
49. NSEL,S,LOC,X,L_HORI
50. F,ALL,FY,10               !力载荷
51. ALLSEL,ALL
52. save

2011-4-11 13:58 上传
下载附件 (39.32 KB)
【提取APDL】
1. !进行一次QRDAMP分析，以生成包含K、M、C和RHS的FULL文件
2. /SOLU
3. ANTYPE,MODAL
4. MODOPT,QRDAMP,2,2
5. SOLVE
6. !将对应的矩阵提取到文件中
7. /AUX2
8. FILE,re,FULL
9. HBMAT,K_RHS,txt,ASCII,,STIFF,YES
10. HBMAT,M,txt,ASCII,,MASS,YES
11. HBMAT,C,txt,ASCII,,DAMP,YES
12. FINISH
【ANSYS谐响应分析】(计算完毕后，手动把受力点的频响结果存在ree.txt中)
1. /SOLU
2. ANTYPE,HARM
3. HARFRQ,0,502
4. NSUBST,251
5. KBC,1   
6. HROPT,FULL  
7. HROUT,OFF   
8. LUMPM,0
9. EQSLV, ,1e-008,
10. SOLVE
【读取&计算】(运行APP_HP_From_ANSYS.py之前在当前目录准备刚才ANSYS计算目录下的K_RHS.txt,M.txt,K.txt,ree.txt)
1. # -*- coding: cp936 -*-
2. #2011.4 重构动力学计算程序
3. #使之可计算线性问题的瞬态和谐响应
4. #使之可从文件读入矩阵
5. #范雨 rainyboy@
6.
7. #本文件包含：一个应用
8. #从ANSYS中导出一个模型的M,K,C和F信息，进行谐响应分析。
9. #与ANSYS的谐响应结果进行对比。
10.
11. import numpy as math;
12. from MechanicMode import *;
13. from HP_FullMethod import *;
14. from ReadMatricFromFile import*;
15.
16. if __name__ == '__main__':
17.     K_RHS_filename = "K_RHS.txt";
18.     M_filename = "M.txt";
19.     C_filename = "C.txt";
20.     df = 2;
21.     freqzone = (0,500);
22.     (K,RHS) = ReadDensMatrixFromFile(K_RHS_filename,True,True);
23.     M = ReadDensMatrixFromFile(M_filename,True,False);
24.     C = ReadDensMatrixFromFile(C_filename,True,False);
25.     #K*1e-5
26.     InitCondition = math.matrix(math.zeros((RHS.shape[1],3),dtype=float));
27.     def forceFunc(f):
28.         return RHS;
29.     Mode = LinearMechanicMode(M,K,C,InitCondition,forceFunc);
30.     Method = FullMethod();
31.     Method.setMode(Mode);
32.     Method.solve(df,freqzone);
33.     #读入ANSYS的求解结果y
34.     f = open("ree.txt");
35.     fcmp = [];
36.     dcomp = [];
37.     while(True):
38.         line = f.readline();
39.         if not line:
40.             break;
41.         infolist = line.split();
42.         fcmp.append(float(infolist[0]));
43.         dcomp.append(float(infolist[1]));
44.     #绘对比图
45.     Method.ComapareFreqHist(fcmp,dcomp,4);
【结果对比图】

2011-4-11 13:58 上传
下载附件 (33.91 KB)

【控制台输出】
正在读取：K_RHS.txt
文件说明：Stiffness matrix from ANSYS FULL file dumped into Harwell-Boeing format         
行指针个数：7
行标个数：12
非零数据个数：12
右端数据个数：6

正在读取：M.txt
文件说明：     Mass matrix from ANSYS FULL file dumped into Harwell-Boeing format         
行指针个数：7
行标个数：12
非零数据个数：12
右端数据个数：6

正在读取：C.txt
文件说明：  Damping matrix from ANSYS FULL file dumped into Harwell-Boeing format         
行指针个数：7
行标个数：12
非零数据个数：12
右端数据个数：6

=============================================================
======5.值得注意的事情============================================
=============================================================
1，为什么在提取矩阵之前要进行QRDAMP的模态分析？
ANSYS帮助说，BHMAT命令中来获取DAMP参数，仅当进行考虑阻尼的模态分析时才有效：
DAMP  —  Write damping matrix to output matrix file. Only valid for damped modal analyses.
因此，我选择了QRDAMP，我只是借助该方法生成FULL文件，并不关注其求解结果

2，注意，在本例中（以及大多数场合），ANSYS导出的矩阵是对称的，即只导出了下三角（或上三角），在ReadFromFile.py中的ReadDensMatrixFromFile函数作了相应的处理，使得到的矩阵是一个对称的满阵。而ReadSparseMatrixFromFile函数则没有。

3，注意，导出的矩阵已经引入了约束条件，即固定的自由度已经被删去了，这就是为什么上述测试用例中，明明是3个BEAM3单元（每个单元3个自由度），得到的矩阵只有6*6。

补充一下，如何得到矩阵的各行与模型各自由度之间的映射关系：
在HBMAT命令中打开mapping开关，即形如（最后那个YES即是此开关）：
HBMAT,KK_RHS,txt,ASCII,,STIFF,YES,YES就会得到名为KK_RHS.mapping的文件，其内容大概为：    Matrix Eqn          Node    DOF     
             1             3    UX  
             2             3    UY  
             3             3    ROTZ
             4             2    UX  
             5             2    UY  
             6             2    ROTZ

通过这些内容就可以知道所关注的自由度对应的是哪个方程（矩阵的哪一行/列）了。

其原理很简单，即使用ansys的超单元即可解决问题。定义超单元，然后列出超单元的刚度矩阵即可。
面是一个小例题，自可明白。
/prep7
k,1
k,2,3000
l,1,2
et,1,beam3
mp,ex,1,2e5
mp,prxy,1,0.3
r,1,5000,2e7,200
lesize,all,,,10
lmesh,all
finish
!----以上正常建立模型，不必施加约束和荷载
/solu
antype,7 !substructuring分析类型
seopt,matname,1 !设置文件名称和刚度矩阵类型(刚度，质量，阻尼等)
nsel,all !选择所有节点
m,all,all !定义所有节点自由度为主自由度
solve !求解
selist,matname,3 !列出整体刚度矩阵

打印 | 推荐 | 订阅 | 收藏
ansys刚度矩阵Harwell-Boeing格式的具体含义讨论
窗体顶端

hefeng902
#1
2008-2-22 11:09   

ANSYS8中有一个新功能，就是提供了一个直接提取刚度矩阵的方法，其存储格式为Harwell-Boeing格式，不知是什么意思，看英文说明看了半天也没看懂，求高手说明一下（最好用中解释一下），谢谢。刚度矩阵如下：

Stiffness matrix from ANSYS FULL file dumped into Harwell-Boeing format         
            25             5             8             8             4
RSA                        4             4             8             0
(I14)           (I14)           (d25.15)            (d25.15)            
F                          1             4
             1
             3
             6
             8
             9
             3
             1
             3
             4
             2
             4
             3
             4
    0.300000000000000D+09
    0.312000000000000D+10
    0.300000000000000D+09
    0.300000000000000D+09
    0.312000000000000D+10
    0.500000000000000D+09
    0.200000000000000D+10
    0.100000000000000D+10
    0.000000000000000D+00
    0.000000000000000D+00
    0.000000000000000D+00
    0.000000000000000D+00

hefeng902 修改于2008-02-22 11:37

iamperisher

#2
2008-2-22 15:29   

HBMAT命令法提取整体矩阵
⑴ 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命令得到质量矩阵文件。
⑵ Harwell-Boeing文件格式
        用HBMAT命令可输出结构刚度矩阵、质量矩阵和阻尼矩阵，其文件记录格式为大型稀疏矩阵的标准交换格式，采用索引存储方法仅记录矩阵的非零元素。文件基本格式是前面有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，并输出右边项。
摘自【ANSYS工程结构数值分析】
书中还详细讲了刚度矩阵的提取与处理等。

hefeng902

积分 166
帖子 114

#3
2008-2-22 16:46   

首先强烈地向iamperisher致敬，回答得很详细，很好，谢谢。对照起来看，受益很大。
操作中还有些问题望你再指点一下，如何将稀疏矩阵转为普通的矩阵呢，是下面这步吗？
采用命令/AUX2
$FILE,mywork,full$HBMAT,mystiff,txt,ASCII,STIFF,YES$FINISH将二进制mywork.full文件输出为ASCII码文件mystiff.txt，并输出右边项。

这几步如何操作啊，好象没有见过前输$的，这一步能不能详细一点呢？

hefeng902

积分 166
帖子 114

#4
2008-2-24 14:29   

通过摸索，大致搞清了以上问题，为了更少了走老路，把摸索的结果大致介绍一下：
iamperisher回答中
$FILE,mywork,full$HBMAT,mystiff,txt,ASCII,STIFF,YES$FINISH将二进制mywork.full文件输出为ASCII码文件mystiff.txt，并输出右边项。
的实现：FILE, Fname, Ext, --

Specifies the data file where results are to be found.

POST1: Set Up

POST26: Set Up

MP ME ST DY <> PR EM <> FL PP ED

Fname
File name and directory path (248 characters maximum, including directory). If you do not specify a directory path, it will default to your working directory and you can use all 248 characters for the file name.

The file name defaults to Jobname.

Ext
Filename extension (8 character maximum).

If Fname is blank, the extension defaults to RST (for structural, fluid, or coupled-field analyses), to RTH (for thermal or electrical analyses), to RMG (for magnetic analyses), or to RFL (for FLOTRAN analyses). For postprocessing reduced structural analyses in POST26, use the RDSP extension for displacements from transient dynamic analyses or the RFRQ extension from harmonic response analyses.

Menu Paths
Main Menu>General Postproc>Data & File Opts
Main Menu>TimeHist Postpro>Settings>File
Utility Menu>File>List>Binary Files
Utility Menu>List>Files>Binary Files
这一步一定要先读取数据，Main Menu>General Postproc>Data & File Opts 不然导不出来。
$HBM