杭州锐达数字技术有限公司
查看: 21880|回复: 68
打印 上一主题 下一主题

[前后处理] 再论ANSYS中的总体矩阵提取(In Python)

? [复制链接]
跳转到指定楼层
楼主
发表于 2011-4-11 14:10 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
============为什么折腾这个文档========
我有一个计算线性动力学方程组的瞬态、谐响应和静力学的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. #从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);
复制代码




【结果对比图】


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


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









Package.rar

9.92 KB, 阅读权限: 20, 下载次数: 52

全部代码

再论ANSYS中的总体矩阵提取.pdf

394.03 KB, 下载次数: 294

pdf

评分

2

查看全部评分

本帖被以下淘专辑推荐:

沙发
发表于 2011-4-12 16:15 | 只看该作者
非常感兴趣,可惜权限不够,可以发到我信箱里么,谢谢了xuexi215@163.com
板凳
?楼主| 发表于 2011-4-12 18:22 | 只看该作者
补充一下,如何得到矩阵的各行与模型各自由度之间的映射关系:
在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



通过这些内容就可以知道所关注的自由度对应的是哪个方程(矩阵的哪一行/列)了。
地板
发表于 2011-4-14 16:39 | 只看该作者
回复 3 # Rainyboy 的帖子

前辈整理的很全面。
请问:在【提取&计算】部分,有没有Matlab相关计算程序可供参考哦。
谢谢。
5
?楼主| 发表于 2011-4-14 17:07 | 只看该作者
回复 4 # mlcc 的帖子

我手里没有,因为我用python的缘故,不过刚才我看了下MATLAB帮助,没发现CSC形式的稀疏矩阵的直接创建函数,不过可以用Sparse函数或者spconvert函数创建,这样的话需要自己写一些代码了。
6
发表于 2011-5-17 15:34 | 只看该作者
不错,比较感兴趣,但是不知道如何运行呢....
能否给像我这样的python盲一些快速安装和运行程序的指导?呵呵
7
?楼主| 发表于 2011-5-17 17:25 | 只看该作者
回复 6 # 心灯 的帖子

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


我那会儿白手起家的时候记录的……
8
发表于 2011-5-17 23:11 | 只看该作者
本帖最后由 心灯 于 2011-5-17 23:14 编辑

好哎,谢谢啦,真没有想到还有pythonxy这样的套装呢,原来的时候我就是对python的各种plugin等的比较困惑,所以一直没有怎么去学....
9
发表于 2011-5-18 09:11 | 只看该作者
本帖最后由 16443 于 2011-5-18 16:10 编辑


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格式的文本文件。

再次修改:当有右边项时,原程序不能输出右边项。修改以后,当有右边项时,输出矩阵的最后一列即为右边项。
10
?楼主| 发表于 2011-5-19 15:59 | 只看该作者
回复 9 # 16443 的帖子

挺好的说明和补充!
看了你的博客,有很多技术细节哈,还在学校里念博士?
11
发表于 2011-5-19 16:20 | 只看该作者
本帖最后由 16443 于 2011-5-19 16:20 编辑

回复 10 # Rainyboy 的帖子

是,正纠结着呢……

你是北航在读还是毕业 ?
12
?楼主| 发表于 2011-5-19 18:57 | 只看该作者
回复 11 # 16443 的帖子

也是马上就要毕业了……同纠结啊
可以说说你的主要方向是什么呢?
我所在的教研室主要研究的是旋转部件的强度、365bet开户推荐_365bet赔率怎么看_365bet怎么不能存和可靠性问题,一般的同学用ANSYS用得多,但是牛一些的师兄都还是有各自的看家程序……拿ANSYS验证验证就行了……
13
发表于 2011-5-23 09:23 | 只看该作者
回复 12 # Rainyboy 的帖子

我毕业还没希望那,没文章到手。
我主要作一些结合面相关,想和商用有限元耦合,目前还没找到办法。
你们可靠性分析也有自己的程序?
14
?楼主| 发表于 2011-5-23 09:53 | 只看该作者
回复 13 # 16443 的帖子

嗯,疲劳寿命计算有一些自己的程序,经典的高低周(和各种修正)的计算方法,也有引入概率化之后的高低周寿命计算方法。不过疲劳寿命这一领域似乎没什么靠谱的理论……主要还是要大量的试验才行吧……
15
发表于 2011-5-24 11:07 | 只看该作者
回复 14 # Rainyboy 的帖子

对,好像没什么理论性的东西,大部分都是基于实验。
现在好像有可靠性方面的国家课题
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表