再论ANSYS中的总体矩阵提取(In Python)

Rainyboy

============为什么折腾这个文档========
我有一个计算线性动力学方程组的瞬态、谐响应和静力学的python程序，现希望开发一个将ANSYS组集好的总体矩阵导入该PYTHON程序中的接口。
该问题可分解为：
[STEP1]  [ANSYS]->[包含矩阵信息的文件]
[STEP2]  [包含矩阵信息的文件]->[python通用数据对象]  
[STEP3]  [python通用数据对象]->[程序特定数据对象]->[进行计算]
因此检索了一些帖子，基本上完成了这项工作，本文是对[STEP1]和[STEP2]的整理，并且利用[STEP3]对结果进行了验证

============主要内容==================
1，了解从ANSYS中提取总体矩阵和载荷向量的方法；
2，了解提取出来的矩阵是怎样表示的；
3，说明在Python中，如何读取这样的矩阵；
4，构造一个简单的算例，说明整个【建模】-【提取】-【读取】过程及其正确性；

======================================
=========站内检索综述====================
======================================
检索词：提取 矩阵
得到21个结果，代表性的帖子有下面这9个：

编号[1]
标题：ansys中怎样提取质量，刚度，阻尼矩阵？
地址：http://forum.vibunion.com/thread-34163-1-1.html
要点：pengweicai给出了一段网上最常见的提取代码,该程序以fortran写成，可以利用.full文件以及一些列约定将ANSYS中的总体矩阵读入FORTRAN中。

编号[2]
标题：如何得知HBMAT命令提取的质量、刚度矩阵对应的自由度？
地址：http://forum.vibunion.com/thread-83591-1-1.html
要点：提出了使用HBMAT命令提取稀疏矩阵时常见的问题：我们如何知道提取出来的信息是怎么储存的呢？

编号[3]
标题：[分享]ANSYS中整体、单元刚度和质量矩阵的提取
地址：http://forum.vibunion.com/thread-2159-1-1.html
要点：在该帖子的7楼，其实已经给出了帖子[2]中问题的解答，即HBMAT中提取出来的矩阵是Harwell-Boeing格式的，并且给出了该格式的细节，可惜是英文的，没引起多少关注。

编号[4]
标题：帮我看看提取的刚度与质量矩阵
地址：http://forum.vibunion.com/thread-81771-1-1.html
要点：这个帖子所示的矩阵并非是使用HBMAT命令提出出来的，而应该是SELIST命令列举出来的未压缩的矩阵，后续楼层的回帖给了大家一个提示，即有可能提取出来的矩阵是引入了边界条件的（即删除了被约束的行和列的）。

编号[5]
标题：提取刚度矩阵的问题
地址：http://forum.vibunion.com/thread-78639-1-1.html
要点：本帖作者的工作是基于单元刚度矩阵的，因此ANSYS中提取的单元刚度矩阵是否处于总体坐标系就成为问题。该问题并非本文内容，但仍值得关注。

编号[6]
标题：提取刚度矩阵丢失节点的问题
地址：http://forum.vibunion.com/thread-78852-1-1.html
要点：帖子[5]作者的又一帖，在这里帖子[5]的问题得到了欧阳中华老师的回答。

编号[7]
标题：提取刚度矩阵的ANSYS操作过程
地址：http://forum.vibunion.com/thread-61253-1-1.html
要点：实际上这就是使用HBMAT从ANSYS中提取总体矩阵的全过程！只是还有一些细节待确定。

编号[8]
标题：提取整体刚度矩阵、质量矩阵及阻尼矩阵的简单方法
地址：http://forum.vibunion.com/thread-8014-1-1.html
要点：给出了利用“不减缩的”子结构方法来得到总体矩阵的方法（这也是网络上常见的代码之一）

编号[9]
标题：质量矩阵、刚度矩阵如何提取？
地址：http://forum.vibunion.com/thread-56104-1-1.html
要点：16443在5楼的回帖中给出了提取刚度矩阵的三种方法

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

百度检索：提取 矩阵
比较好的帖子有：

编号[10]
来源:百度文库
标题:怎样从ansys中提取单元刚度矩阵与质量矩阵
地址:http://wenku.baidu.com/view/3cf5e567f5335a8102d220d9.html
要点:这应该就是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格式数据的说明
地址:http://forum.simwe.com/archiver/tid-924778.html
要点:比[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

复制代码

【提取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@188.com
   
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);

复制代码

【结果对比图】


【控制台输出】
正在读取：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。

====================================================
欢迎你来算法区讨论python：
http://forum.vibunion.com/forum-25-1.html


附件包含全部代码、可执行模块和测试用例：

leemu215

非常感兴趣，可惜权限不够，可以发到我信箱里么，谢谢了xuexi215@163.com

Rainyboy

补充一下，如何得到矩阵的各行与模型各自由度之间的映射关系：
在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



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

mlcc

回复 3 # Rainyboy 的帖子

前辈整理的很全面。
请问：在【提取&计算】部分，有没有Matlab相关计算程序可供参考哦。
谢谢。

Rainyboy

回复 4 # mlcc 的帖子

我手里没有，因为我用python的缘故，不过刚才我看了下MATLAB帮助，没发现CSC形式的稀疏矩阵的直接创建函数，不过可以用Sparse函数或者spconvert函数创建，这样的话需要自己写一些代码了。

心灯

不错，比较感兴趣，但是不知道如何运行呢....
能否给像我这样的python盲一些快速安装和运行程序的指导？呵呵

Rainyboy

回复 6 # 心灯 的帖子

哦，这样啊，看看这个帖子呗：
Python：一场从零开始的奇妙旅程
http://forum.vibunion.com/forum- ... fromuid-159019.html


我那会儿白手起家的时候记录的……

心灯

好哎，谢谢啦，真没有想到还有pythonxy这样的套装呢，原来的时候我就是对python的各种plugin等的比较困惑，所以一直没有怎么去学....

16443

mlcc 发表于 2011-4-14 16:39

                               
登录/注册后可看大图

回复 3 # Rainyboy 的帖子

前辈整理的很全面。

http://forum.simwe.com/viewthread.php?tid=912201

知道了Harwell-Boeing的矩阵的储存格式以后，就可以将按照Matlab稀疏矩阵的格式转换为稀疏矩阵。
可以采用下面网址的m文件进行上述转换，http://people.sc.fsu.edu/~burkardt/m_src/hb_to_mm/hb_to_mm.m。
但这个m文件只是输出一个文本文件。
为了可以在Matlab中直接应用，对其进行了稍微修改，修改后的命令有3种形式：

kkk=hb2mm('kkk_HB.txt','kkk_mm.txt');
kkk=hb2mm('kkk_HB.txt');
hb2mm('kkk_HB.txt','kkk_mm.txt');

其中，kkk即为Matlab格式的稀疏矩阵。kkk_mm.txt为输出的含有稀疏矩阵的文本文件。kkk_HB.txt为Harwell-Boeing格式的文本文件。

再次修改：当有右边项时，原程序不能输出右边项。修改以后，当有右边项时，输出矩阵的最后一列即为右边项。

Rainyboy

回复 9 # 16443 的帖子

挺好的说明和补充！
看了你的博客，有很多技术细节哈，还在学校里念博士？

16443

回复 10 # Rainyboy 的帖子

是，正纠结着呢……

你是北航在读还是毕业 ？

Rainyboy

回复 11 # 16443 的帖子

也是马上就要毕业了……同纠结啊
可以说说你的主要方向是什么呢？
我所在的教研室主要研究的是旋转部件的强度、振动和可靠性问题，一般的同学用ANSYS用得多，但是牛一些的师兄都还是有各自的看家程序……拿ANSYS验证验证就行了……

16443

回复 12 # Rainyboy 的帖子

我毕业还没希望那，没文章到手。
我主要作一些结合面相关，想和商用有限元耦合，目前还没找到办法。
你们可靠性分析也有自己的程序？

Rainyboy

回复 13 # 16443 的帖子

嗯，疲劳寿命计算有一些自己的程序，经典的高低周（和各种修正）的计算方法，也有引入概率化之后的高低周寿命计算方法。不过疲劳寿命这一领域似乎没什么靠谱的理论……主要还是要大量的试验才行吧……

16443

回复 14 # Rainyboy 的帖子

对，好像没什么理论性的东西，大部分都是基于实验。
现在好像有可靠性方面的国家课题