提取刚度矩阵

============为什么折腾这个文档========我有一个计算线性动力学方程组的瞬态、谐响应和静力学的python程序，现希望开发一个将ANSYS组集好的总体矩阵导入该PYTHON程序中的接口。该问题可分解为：[STEP1][ANSYS]-[包含矩阵信息的文件][STEP2][包含矩阵信息的文件]-[python通用数据对象][STEP3][python通用数据对象]-[程序特定数据对象]-[进行计算]因此检索了一些帖子，基本上完成了这项工作，本文是对[STEP1]和[STEP2]的整理，并且利用[STEP3]对结果进行了验证============主要内容==================1，了解从ANSYS中提取总体矩阵和载荷向量的方法；2，了解提取出来的矩阵是怎样表示的；3，说明在Python中，如何读取这样的矩阵；4，构造一个简单的算例，说明整个【建模】-【提取】-【读取】过程及其正确性；=========站内检索综述====================检索词：提取矩阵得到21个结果，代表性的帖子有下面这9个：编号[1]标题：ansys中怎样提取质量，刚度，阻尼矩阵？地址：要点：pengweicai给出了一段网上最常见的提取代码,该程序以fortran写成，可以利用.full文件以及一些列约定将ANSYS中的总体矩阵读入FORTRAN中。编号[2]标题：如何得知HBMAT命令提取的质量、刚度矩阵对应的自由度？地址：要点：提出了使用HBMAT命令提取稀疏矩阵时常见的问题：我们如何知道提取出来的信息是怎么储存的呢？编号[3]标题：[分享]ANSYS中整体、单元刚度和质量矩阵的提取地址：要点：在该帖子的7楼，其实已经给出了帖子[2]中问题的解答，即HBMAT中提取出来的矩阵是Harwell-Boeing格式的，并且给出了该格式的细节，可惜是英文的，没引起多少关注。编号[4]标题：帮我看看提取的刚度与质量矩阵地址：要点：这个帖子所示的矩阵并非是使用HBMAT命令提出出来的，而应该是SELIST命令列举出来的未压缩的矩阵，后续楼层的回帖给了大家一个提示，即有可能提取出来的矩阵是引入了边界条件的（即删除了被约束的行和列的）。编号[5]标题：提取刚度矩阵的问题地址：要点：本帖作者的工作是基于单元刚度矩阵的，因此ANSYS中提取的单元刚度矩阵是否处于总体坐标系就成为问题。该问题并非本文内容，但仍值得关注。编号[6]标题：提取刚度矩阵丢失节点的问题地址：要点：帖子[5]作者的又一帖，在这里帖子[5]的问题得到了欧阳中华老师的回答。编号[7]标题：提取刚度矩阵的ANSYS操作过程地址：要点：实际上这就是使用HBMAT从ANSYS中提取总体矩阵的全过程！只是还有一些细节待确定。编号[8]标题：提取整体刚度矩阵、质量矩阵及阻尼矩阵的简单方法地址：要点：给出了利用“不减缩的”子结构方法来得到总体矩阵的方法（这也是网络上常见的代码之一）编号[9]标题：质量矩阵、刚度矩阵如何提取？地址：要点：16443在5楼的回帖中给出了提取刚度矩阵的三种方法=======站外检索略述========================百度检索：提取矩阵比较好的帖子有：编号[10]来源:百度文库标题:怎样从ansys中提取单元刚度矩阵与质量矩阵地址:要点:这应该就是16443在帖子[9]中回复的内容了，全面的总结了在帖子[3,4,5,9]中涉及的问题。编号[11]来源:中华钢结构标题:ansys刚度矩阵Harwell-Boeing格式的具体含义讨论地址:=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格式，也叫CompressedSparseColumn格式。其具体压缩和还原方式见帖子[3]（English）或[11][12](中文)=======3.如何在Python中读入BH格式的矩阵===============================上文说过，Harwell-boeing格式，也叫CompressedSparseColumn格式，而Python.scipy中就有这样的稀疏矩阵：1.classscipy.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.RainyboyTestingMatrixinBHformat2.***46603.RRA******04.(I14)(I14)(d25.15)(d25.15)5.06.27.38.69.010.211.212.013.114.215.116.217.318.419.520.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.FINISH2./CLEAR3./TITLE,CASESTUDY_BEAM_BEAM3BYRAINYBOY4./PREP75./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,038.L,1,2!建立几何体39.ALLSEL,ALL40.LESIZE,ALL,,,MESHCOUNT!设置线段分网数41.LMESH,ALL!分网42.!**********************43.!几何约束44.!**********************45.ALLSEL,ALL46.NSEL,S,LOC,X,0!选择固定端节点47.D,ALL,ALL!设置为约束所有自由度48.ALLSEL,ALL49.NSEL,S,LOC,X,L_HORI50.F,ALL,FY,10!力载荷51.ALLSEL,ALL52.save【提取APDL】1.!进行一次QRDAMP分析，以生成包含K、M、C和RHS的FULL文件2./SOLU3.ANTYPE,MODAL4.MODOPT,QRDAMP,2,25.SOLVE6.!将对应的矩阵提取到文件中7./AUX28.FILE,re,FULL9.HBMAT,K_RHS,txt,ASCII,,STIFF,YES10.HBMAT,M,txt,ASCII,,MASS,YES11.HBMAT,C,txt,ASCII,,DAMP,YES12.FINISH【ANSYS谐响应分析】(计算完毕后，手动把受力点的频响结果存在ree.txt中)1./SOLU2.ANTYPE,HARM3.HARFRQ,0,5024.NSUBST,2515.KBC,16.HROPT,FULL7.HROUT,OFF8.LUMPM,09.EQSLV,,1e-008,10.SOLVE【读