声振论坛 首页 工程软件 查看内容

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

2011-4-11 16:26| 发布者: 雪缘| 查看: 21820| 评论: 0|来自: 振动论坛

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

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

======================================
=========站内检索综述====================
======================================

检索词:提取 矩阵
得到21个结果,代表性的帖子有下面这9个:

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

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

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

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

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

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

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

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

编号[9]
标题:质量矩阵、刚度矩阵如何提取?
地址:http://www.chinavib.com/forum-vi ... fromuid-159019.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
复制代码



[attach]50887[/attach]



【提取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. #从ANSYS中导出一个模型的M,K,C和F信息,进行谐响应分析。
  8. #与ANSYS的谐响应结果进行对比。

  9. import numpy as math;
  10. from MechanicMode import *;
  11. from HP_FullMethod import *;
  12. from ReadMatricFromFile import*;

  13. if __name__ == '__main__':
  14.     K_RHS_filename = "K_RHS.txt";
  15.     M_filename = "M.txt";
  16.     C_filename = "C.txt";
  17.     df = 2;
  18.     freqzone = (0,500);
  19.     (K,RHS) = ReadDensMatrixFromFile(K_RHS_filename,True,True);
  20.     M = ReadDensMatrixFromFile(M_filename,True,False);
  21.     C = ReadDensMatrixFromFile(C_filename,True,False);
  22.     #K*1e-5
  23.     InitCondition = math.matrix(math.zeros((RHS.shape[1],3),dtype=float));
  24.     def forceFunc(f):
  25.         return RHS;
  26.     Mode = LinearMechanicMode(M,K,C,InitCondition,forceFunc);
  27.     Method = FullMethod();
  28.     Method.setMode(Mode);
  29.     Method.solve(df,freqzone);
  30.     #读入ANSYS的求解结果y
  31.     f = open("ree.txt");
  32.     fcmp = [];
  33.     dcomp = [];
  34.     while(True):
  35.         line = f.readline();
  36.         if not line:
  37.             break;
  38.         infolist = line.split();
  39.         fcmp.append(float(infolist[0]));
  40.         dcomp.append(float(infolist[1]));
  41.     #绘对比图
  42.     Method.ComapareFreqHist(fcmp,dcomp,4);
复制代码




【结果对比图】
[attach]50888[/attach]

【控制台输出】
正在读取: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://www.chinavib.com/forum-forumdisplay-fid-25-filter-typeid-typeid-41.html





附件包含全部代码、可执行模块和测试用例:
[attach]50889[/attach]








本文内容由 Rainyboy 提供



[attach]50890[/attach]

[attach]50891[/attach]

最新评论

返回顶部