再论ANSYS中的总体矩阵提取(In Python)
============为什么折腾这个文档========我有一个计算线性动力学方程组的瞬态、谐响应和静力学的python程序，现希望开发一个将ANSYS组集好的总体矩阵导入该PYTHON程序中的接口。
该问题可分解为：
->[包含矩阵信息的文件]
[包含矩阵信息的文件]->
->[程序特定数据对象]->[进行计算]
因此检索了一些帖子，基本上完成了这项工作，本文是对和的整理，并且利用对结果进行了验证
============主要内容==================
1，了解从ANSYS中提取总体矩阵和载荷向量的方法；
2，了解提取出来的矩阵是怎样表示的；
3，说明在Python中，如何读取这样的矩阵；
4，构造一个简单的算例，说明整个【建模】-【提取】-【读取】过程及其正确性；
======================================
=========站内检索综述====================
======================================
检索词：提取 矩阵
得到21个结果，代表性的帖子有下面这9个：
编号
标题：ansys中怎样提取质量，刚度，阻尼矩阵？
地址：http://forum.vibunion.com/thread-34163-1-1.html
要点：pengweicai给出了一段网上最常见的提取代码,该程序以fortran写成，可以利用.full文件以及一些列约定将ANSYS中的总体矩阵读入FORTRAN中。
编号
标题：如何得知HBMAT命令提取的质量、刚度矩阵对应的自由度？
地址：http://forum.vibunion.com/thread-83591-1-1.html
要点：提出了使用HBMAT命令提取稀疏矩阵时常见的问题：我们如何知道提取出来的信息是怎么储存的呢？
编号
标题：[分享]ANSYS中整体、单元刚度和质量矩阵的提取
地址：http://forum.vibunion.com/thread-2159-1-1.html
要点：在该帖子的7楼，其实已经给出了帖子中问题的解答，即HBMAT中提取出来的矩阵是Harwell-Boeing格式的，并且给出了该格式的细节，可惜是英文的，没引起多少关注。
编号
标题：帮我看看提取的刚度与质量矩阵
地址：http://forum.vibunion.com/thread-81771-1-1.html
要点：这个帖子所示的矩阵并非是使用HBMAT命令提出出来的，而应该是SELIST命令列举出来的未压缩的矩阵，后续楼层的回帖给了大家一个提示，即有可能提取出来的矩阵是引入了边界条件的（即删除了被约束的行和列的）。
编号
标题：提取刚度矩阵的问题
地址：http://forum.vibunion.com/thread-78639-1-1.html
要点：本帖作者的工作是基于单元刚度矩阵的，因此ANSYS中提取的单元刚度矩阵是否处于总体坐标系就成为问题。该问题并非本文内容，但仍值得关注。
编号
标题：提取刚度矩阵丢失节点的问题
地址：http://forum.vibunion.com/thread-78852-1-1.html
要点：帖子作者的又一帖，在这里帖子的问题得到了欧阳中华老师的回答。
编号
标题：提取刚度矩阵的ANSYS操作过程
地址：http://forum.vibunion.com/thread-61253-1-1.html
要点：实际上这就是使用HBMAT从ANSYS中提取总体矩阵的全过程！只是还有一些细节待确定。
编号
标题：提取整体刚度矩阵、质量矩阵及阻尼矩阵的简单方法
地址：http://forum.vibunion.com/thread-8014-1-1.html
要点：给出了利用“不减缩的”子结构方法来得到总体矩阵的方法（这也是网络上常见的代码之一）
编号
标题：质量矩阵、刚度矩阵如何提取？
地址：http://forum.vibunion.com/thread-56104-1-1.html
要点：16443在5楼的回帖中给出了提取刚度矩阵的三种方法
========================================
=======站外检索略述========================
========================================
百度检索：提取 矩阵
比较好的帖子有：
编号
来源:百度文库
标题:怎样从ansys中提取单元刚度矩阵与质量矩阵
地址:http://wenku.baidu.com/view/3cf5e567f5335a8102d220d9.html
要点:这应该就是16443在帖子中回复的内容了，全面的总结了在帖子中涉及的问题。
编号
来源:中华钢结构
标题:ansys刚度矩阵Harwell-Boeing格式的具体含义讨论
地址:http://okok.org/forum/viewthread.php?tid=184007
要点:如题，后续楼层给出了一些将矩阵读入ANSYS的APDL（好不容易读出来，又读进去干嘛呢……）
编号
来源:simwe
标题:关于ANSYS（质量、刚度、阻尼）矩阵Harwell-boeing格式数据的说明
地址:http://forum.simwe.com/archiver/tid-924778.html
要点:比更透彻的HB格式说明!
=============================================================
=======1.从ANSYS中提取总体矩阵的方法=================================
=============================================================
1,用/DEBUG命令
2,子结构法
3,HBMAT
详见帖子
PS.个人感觉HBMAT方法最靠谱，一是它的格式（Harwell-boeing）在很多场合都是通用的，二是BHMAT命令是文档化的、功能就是用来提取总体刚度矩阵的命令。因此，相比于子结构法的剑走偏锋，/DEBUG命令的繁复，HBMAT命令方法更“标准”一些，因此在后文只关注此方法。
=============================================================
=======2.BH格式的矩阵是如何表示的===================================
=============================================================
HBMAT命令并不是很复杂的命令，稍复杂的地方是采用该命令提取出来的矩阵是经过压缩的，称为Harwell-boeing格式，也叫Compressed Sparse Column格式。
其具体压缩和还原方式见帖子（English）或(中文)
=============================================================
=======3.如何在Python中读入BH格式的矩阵===============================
=============================================================
上文说过，Harwell-boeing格式，也叫Compressed Sparse Column格式，而Python.scipy中就有这样的稀疏矩阵：
class scipy.sparse.csc_matrix(arg1, shape=None, dtype=None, copy=False, dims=None, nzmax=None)
可以通过HB文件中直接读取的行标指针，行标和数据创建，例如：>>> indptr = array()
>>> indices = array()
>>> data = array()
>>> csc_matrix( (data,indices,indptr), shape=(3,3) ).todense()
matrix([,
,
])
对应的HB文件应为(*号部分表示并非本例关注的数据)：
Rainyboy Testing Matrix in BH format
*** 4 6 6 0
RRA ** ** ** 0
(I14) (I14) (d25.15) (d25.15)
0
2
3
6
0
2
2
0
1
2
1
2
3
4
5
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】
FINISH
/CLEAR
/TITLE,CASE STUDY _BEAM _BEAM3 BY RAINYBOY
/PREP7
/ESHAPE,1 !显示壳单元厚度
!**********************
!几何参数表
!**********************
*SET,L_HORI,0.1 !横梁的长度
*SET,TA,0.005 !正方形截面的边长
*SET,MESHCOUNT,2 !每段的分网数
*SET,IZZ,TA*TA*TA*TA/12 !转动惯量
*SET,IYY,TA*TA*TA*TA/12 !转动惯量
!**********************
!材料参数表
!**********************
*SET,MEX,1.78E11 !弹性模量
*SET,MPRXY,0.3 !泊松比
*SET,MDENS,7850 !密度
!**********************
!相关设置
!**********************
MP,EX,1,MEX !设置材料弹性模量
MP,PRXY,1,MPRXY !设置材料泊松比
MP,DENS,1,MDENS !设置材料密度
BETAD,1E-5 !BETA阻尼
ET,1,BEAM3 !设置平面梁单元
R,1,TA*TA,IZZ,TA !设置截面参数
!DMPRAT,0.10000 !阻尼比
!**********************
!几何->分网
!**********************
TYPE,1 !指定分网类型
MAT,1 !指定材料类型
REAL,1 !指定实参数
K,1,0,0,0 !建立三个关键点
K,2,L_HORI,0,0
L,1,2 !建立几何体
ALLSEL,ALL
LESIZE,ALL,,,MESHCOUNT !设置线段分网数
LMESH,ALL !分网
!**********************
!几何约束
!**********************
ALLSEL,ALL
NSEL,S,LOC,X,0 !选择固定端节点
D,ALL,ALL !设置为约束所有自由度
ALLSEL,ALL
NSEL,S,LOC,X,L_HORI
F,ALL,FY,10 !力载荷
ALLSEL,ALL
save
【提取APDL】
!进行一次QRDAMP分析，以生成包含K、M、C和RHS的FULL文件
/SOLU
ANTYPE,MODAL
MODOPT,QRDAMP,2,2
SOLVE
!将对应的矩阵提取到文件中
/AUX2
FILE,re,FULL
HBMAT,K_RHS,txt,ASCII,,STIFF,YES
HBMAT,M,txt,ASCII,,MASS,YES
HBMAT,C,txt,ASCII,,DAMP,YES
FINISH
【ANSYS谐响应分析】(计算完毕后，手动把受力点的频响结果存在ree.txt中)
/SOLU
ANTYPE,HARM
HARFRQ,0,502
NSUBST,251
KBC,1
HROPT,FULL
HROUT,OFF
LUMPM,0
EQSLV, ,1e-008,
SOLVE
【读取&计算】(运行APP_HP_From_ANSYS.py之前在当前目录准备刚才ANSYS计算目录下的K_RHS.txt,M.txt,K.txt,ree.txt)
# -*- coding: cp936 -*-
#2011.4 重构动力学计算程序
#使之可计算线性问题的瞬态和谐响应
#使之可从文件读入矩阵
#范雨 rainyboy@188.com
#本文件包含：一个应用
#从ANSYS中导出一个模型的M,K,C和F信息，进行谐响应分析。
#与ANSYS的谐响应结果进行对比。
import numpy as math;
from MechanicMode import *;
from HP_FullMethod import *;
from ReadMatricFromFile import*;
if __name__ == '__main__':
K_RHS_filename = "K_RHS.txt";
M_filename = "M.txt";
C_filename = "C.txt";
df = 2;
freqzone = (0,500);
(K,RHS) = ReadDensMatrixFromFile(K_RHS_filename,True,True);
M = ReadDensMatrixFromFile(M_filename,True,False);
C = ReadDensMatrixFromFile(C_filename,True,False);
#K*1e-5
InitCondition = math.matrix(math.zeros((RHS.shape,3),dtype=float));
def forceFunc(f):
return RHS;
Mode = LinearMechanicMode(M,K,C,InitCondition,forceFunc);
Method = FullMethod();
Method.setMode(Mode);
Method.solve(df,freqzone);
#读入ANSYS的求解结果y
f = open("ree.txt");
fcmp = [];
dcomp = [];
while(True):
line = f.readline();
if not line:
break;
infolist = line.split();
fcmp.append(float(infolist));
dcomp.append(float(infolist));
#绘对比图
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
附件包含全部代码、可执行模块和测试用例：
非常感兴趣，可惜权限不够，可以发到我信箱里么，谢谢了xuexi215@163.com 补充一下，如何得到矩阵的各行与模型各自由度之间的映射关系：
在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
通过这些内容就可以知道所关注的自由度对应的是哪个方程（矩阵的哪一行/列）了。
回复 3 # Rainyboy 的帖子
前辈整理的很全面。{:{39}:}
请问：在【提取&计算】部分，有没有Matlab相关计算程序可供参考哦。
谢谢。 回复 4 # mlcc 的帖子
我手里没有，因为我用python的缘故，不过刚才我看了下MATLAB帮助，没发现CSC形式的稀疏矩阵的直接创建函数，不过可以用Sparse函数或者spconvert函数创建，这样的话需要自己写一些代码了。 不错，比较感兴趣，但是不知道如何运行呢....
能否给像我这样的python盲一些快速安装和运行程序的指导？呵呵 回复 6 # 心灯 的帖子
哦，这样啊，看看这个帖子呗：
Python：一场从零开始的奇妙旅程
http://forum.vibunion.com/forum-viewthread-tid-96961-fromuid-159019.html
我那会儿白手起家的时候记录的…… 本帖最后由 心灯 于 2011-5-17 23:14 编辑
好哎，谢谢啦，真没有想到还有pythonxy这样的套装呢，原来的时候我就是对python的各种plugin等的比较困惑，所以一直没有怎么去学....{:{39}:} 本帖最后由 16443 于 2011-5-18 16:10 编辑
mlcc 发表于 2011-4-14 16:39 static/image/common/back.gif
回复 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格式的文本文件。
再次修改：当有右边项时，原程序不能输出右边项。修改以后，当有右边项时，输出矩阵的最后一列即为右边项。 回复 9 # 16443 的帖子
挺好的说明和补充！
看了你的博客，有很多技术细节哈，还在学校里念博士？ 本帖最后由 16443 于 2011-5-19 16:20 编辑
回复 10 # Rainyboy 的帖子
是，正纠结着呢……
你是北航在读还是毕业 ？ 回复 11 # 16443 的帖子
也是马上就要毕业了……同纠结啊
可以说说你的主要方向是什么呢？
我所在的教研室主要研究的是旋转部件的强度、振动和可靠性问题，一般的同学用ANSYS用得多，但是牛一些的师兄都还是有各自的看家程序……拿ANSYS验证验证就行了…… 回复 12 # Rainyboy 的帖子
我毕业还没希望那，没文章到手。
我主要作一些结合面相关，想和商用有限元耦合，目前还没找到办法。
你们可靠性分析也有自己的程序？ 回复 13 # 16443 的帖子
嗯，疲劳寿命计算有一些自己的程序，经典的高低周（和各种修正）的计算方法，也有引入概率化之后的高低周寿命计算方法。不过疲劳寿命这一领域似乎没什么靠谱的理论……主要还是要大量的试验才行吧…… 回复 14 # Rainyboy 的帖子
对，好像没什么理论性的东西，大部分都是基于实验。
现在好像有可靠性方面的国家课题