Python-Matlab联合编程Abaqus高级后处理技术研究

任海峰，高　鸣

(海军航空工程学院，山东 烟台　264001)

Python-Matlab联合编程Abaqus高级后处理技术研究

任海峰，高鸣

(海军航空工程学院，山东 烟台264001)

摘要：有限元技术的深度应用对通用软件计算结果输出方式的灵活性、定制性提出了更高的要求；针对大型通用有限元软件Abaqus批量输出节点应力较为困难的问题，提出利用Python-Matlab编程定制输出Abaqus节点应力的新方法，分别用Abaqus应力重构计算输出法和Abaqus节点应力内核调用输出法实现了固体火箭发动健康监测所使用的粘接应力传感器节点应力计算值的批量定制输出；结果表明：两种方法正确有效，精度较高，灵活性好，解释性好。

关键词：固体火箭发动机；粘接应力传感器；Python-Matlab编程；批量输出; Abaqus高级后处理；重构计算; 内核调用

有限元方法作为高精度的计算技术在工程中得到了广泛使用，目前Abaqus、Ansys等通用有限元软件已经成为数值分析的有力工具。但是对有限元结果的深入利用和对软件的二次开发都涉及对计算结果更为灵活的访问和使用。而这恰恰是大型封装通用有限元软件不足。以Abaqus为例，单元应力以分量形式储存在单元积分点，而非节点[1-3]。在研究中常常需要批量输出节点应力和坐标，以便进行仿真结果的对比分析。此外，一些仿真过程中由于模型的改变，重新划分网格，导致网格节点和单元的重新排列，网格节点和单元的数目、位置不同，这给仿真结果的分析带来极大不便。然而Abaqus提供的节点应力输出方式不适合批量输出，需要寻求新的解决方案。但由于软件自身的封闭性，多平台交互的困难性和节点应力计算输出的复杂性，尚未见到介绍Abaqus节点应力批量输出的文献。

Abaqus提供了Python脚本接口，利用被誉为“胶水语言”的Python语言的可扩充性、可移植性、可嵌入性、可解释性，通过Python编程实现Abaqus二次开发，可以实现Abaqus/CAE数据交互难于实现的问题，例如参数化建模、自动化分析、智能化后处理和批量化输入输出等[2,4]。

本文通过原理分析、变化实现、实例对比为研究者提供基于Python和Abaqus高级后处理的节点应力批量输出解决方案，为软件的深度开发利用提供了便捷的技术途径。

1问题描述

采用Python编程定制输出节点应力有两种方法：Abaqus节点应力重构计算输出法和Abaqus节点应力内核调用输出法。前者利用Python编程提取几何集单元应力分量、单元、节点信息，利用Excel处理输出数据格式，利用Matlab实现由单元应力向节点应力的插值外推。后者利用Python编程自动生成节点列表，并调用Abaqus节点输出内核，用Excel处理输出数据格式，用Matlab实现后处理。

为便于分析对比，本文以固体火箭发动机上使用的微型粘接应力传感器分析结果的后处理为例进行对比分析。如图1所示，某E型应力传感器膜盒为例，膜盒零件Part-1，膜盒底面为几何集Set-1，(在AbaqusCAEPart模块，Tools菜单→Set→CreateSet定义膜盒底面为Set-1(type=Geometry))。仿真结果odb文件为Job-31.Odb，分析步为Step-1，单元类型为C3D8R，要求输出Set-1所有节点应力。

图1　E型传感器膜盒及几何集Set-1示意图

2Abaqus节点应力重构计算输出法

根据有限元的近似解性质,应力和应变近似解一定在精确解上下震荡的，但某些点上的解正好和精确解相等，即最佳应力点[5]。根据以应力为自变量的最小位能原理，采用高斯数值积分，由高斯积分的性质，可知积分点处为应力的精确解。由泰勒级数展开式可知，数值计算中积分比微分计算精度高[1,6]，因此，Abaqus中的应力存储在高斯积分点，而不是节点，节点应力由单元应力插值外推[1]。因此，要实现Abaqus高级后处理批量输出节点应力，首先要输出几何集单元应力分量、单元、节点信息，然后采用计算单元应力值，如最大主应力，等效应力(Mises应力)，再搜索节点相关单元，利用插值外推求解节点应力。其间不同软件数据格式的要求不同，需要进行数据格式的转换。

2.1　Abaqus数据结构

为了编程实现节点应力的自动输出必须全面掌握Abaqus的数据结构，以便实现输出过程。Abaqus主要有SESSION，ODB，MDB三种对象[1,7-10]，本方法使用ODB，MDB二种对象，其数据结构如下：

2.1.1AbaqusODB对象数据结构

如图2所示，ODB对象主要对包含计算模型对象数据(ModelData)和计算结果数据(ResultData)，计算模型对象数据包含装配体(rootAssenbly)、零件(parts)、界面分类(sectionCategories)、材料 (materials)等子对象，计算结果数据steps包含分析步(step)、帧(frame)、历史变量输出(historyoutputs)和场变量输出(fieldoutputs)等。

场变量的读取路径：odb.setps[].frames[].fieldOutputs[]；

场变量包括：应力分量—′S′；应变—′E′；位移—′U′；

历史变量的读取路径：odb.setps[].historyRegions[].historyOutputs[]

2.1.2AbaqusMDB对象数据结构

如图3所示，MDB对象主要对包含计算模型对象数据(ModelData)和工作对象数据(JobData)，其中model对象，包含零件(parts)、材料(materials)等子对象。零件(parts)子对象包含单元(elements)、集合(sets)和节点(nodes)子对象，其下又包含编号(label)、坐标(coordinate)、单元所属节点组(connectivity)等属性。

几何集的读取路径：

mdb.models[].parts[].sets[]

odb.rootAssembly.instances[].elementSets[]

几何集所属单元、节点信息读取路径：

单元编号：

mdb.models[].parts[].sets[].elements[].elementLabel

odb.rootAssembly.instances[].elementSets[].elements[].elementLabel

单元所属节点组：

mdb.models[].parts[].sets[].elements[].connectivity

odb.rootAssembly.instances[].elementSets[].elements[].connectivity

节点坐标：

mdb.models[].parts[].sets[].nodes[].coordinates

图2　ODB对象数据结构

图3　MDB对象数据结构

2.2　Abaqus节点应力重构计算输出法程序源码及说明

采用Python编程以读取场变量的方式，在Abaqus计算结果文件Job-31.odb中读取Set-1所有单元应力分量并写入Scomponent.txt。读取Set-1的每个单元所属节点编号写入elmementofset.txt。读取Set-1的节点坐标写入nodeofset.txt。程序源码及说明如下:

算法1几何集单元应力分量、单元信息和节点信息输出算法

#导入odbAccess模块

fromodbAccessimport*

#打开odb文件

myodb=openOdb(′Job-31.odb′)

#读取Step-1中frames[-1]中应力场

stressfield=myodb.steps[′Step-1′].frames[-1].fieldOutputs[′S′]

#定义几何集Set-1

skinset=myodb.rootAssembly.instances[′PART-1-1′].elementSets[′SET-1′]

#取几何集Set-1的应力场

field1=stressfield.getSubset(region=skinset,position=INTEGRATION_POINT,elementType= ′C3D8R′)

#读取几何集Set-1的所有单元的应力分量数据

val=field1.values

#打开Scomponent.txt的并将Set-1所有单元编号及应力分量数据按行循环写入文件，完成后屏幕输出ok，并关闭Scomponent.txt文件

withopen(′D:/Scomponent.txt′,′w′)asf1:

forvinval:

f1.writelines(str(v.elementLabel)+′ ′)

f1.writelines(str(v.data)+′ ′)

else:

print′ok′

f1.close

#读取几何集每个单元所属节点组节点编号

#打开elmementofset.txt的并将Set-1所有单元组成节点编号按行循环#写入文件，完成后屏幕输出ok2，并关闭elmementofset.txt文件

withopen(′D:/elmementofset.txt′,′w′)asf2:

forvinval:

elementnode=myodb.rootAssembly.instances[′PART-1-1′].elements[v.elementLabel-1].connectivity

f2.writelines(str(v.elementLabel)+′ ′)

f2.writelines(str(elementnode)+′ ′)

else:

f2.close

print′ok2′

#读取Set-1所属节点坐标，并将节点坐标按行写入文本文件#nodeofset.txt，完成后屏幕输出ok3，并关闭nodeofset.txt文件

nd=mdb.models[′Model-1′].parts[′Part-1′].sets[′Set-1′].nodes

withopen(′D:/nodeofset.txt′,′w′)asf3:

v=1

whilev<=len(nd):

nodescoors=mdb.models[′Model-1′].parts[′Part-1′].

sets[′Set-1′].nodes[v-1].coordinates

nodeslabel=mdb.models[′Model-1′].parts[′Part-1′].

sets[′Set-1′].nodes[v-1].label

f3.writelines(str(nodeslabel)+′′)

f3.writelines(str(nodescoors)+′ ′)

v=v+1

else:

f3.close

print′ok3′

#关闭odb文件

myodb.close()

2.3　数据格式处理

Python编程输出的数据虽然被保存为文本文件，但其保持了Abaqus原数据格式，含有逗号，方括号，圆括号等符号，为了方便Matlab读入处理，需要用Excel进行数据格式的整理转化。数据转换方法：采用Excel打开以上Python输出文件*.txt文件。采用Excel文本导入向导功能删除数据符号。采用Excel数据分列向导功能删除剩余符号；保存数据为*.xls，*.txt等Matlab可载入的数据文件。

2.4　单元应力重构计算

Abaqus存储为单元积分点的应力分量,例如，三维问题的应力分量为S11，S22，S33，S12，S13，S23，其中S11，S22，S33分别为xx,yy,zz方向正应力，S12为xy方向剪应力，S13为xz方向剪应力，S23为yz方向剪应力，由于xy,xz,yz与yx,zx,zy方向两两对应相等[1,6]，则单元应力矩阵为

(1)

单元主应力为单元应力矩阵的特征值σ1,σ2,σ3；单元等效应力(VonMises应力)[6]：

(2)

算法2单元应力重构算法

clear;clc%清屏

loadD: eport.txt；%载入应力分量数据

fori=1:1:length(Scomponent);%读入单元应力分量

s11=Scomponent(i,2);

s22=Scomponent(i,3);

s33=Scomponent(i,4);

s12=Scomponent(i,5);

s13=Scomponent(i,6);

s23=Scomponent(i,7);

A=[s11s12s13;s12s22s23;s13s23s33]; %构造单元应力矩阵

ps=eig(A);%求解主应力

mpofset(i,1)=Scomponent(i,1); %写入单元编号到mpofset

mpofset(i,2)=ps(1); %计算写入单元主应力到mpofset

mpofset(i,3)=ps(2);

mpofset(i,4)=ps(3);

mpofset(i,5)=max(ps); %写入单元最大主应力到mpofset

mpofset(i,6)=sqrt(0.5*((ps(1)-ps(2))^2+(ps(2)-ps(3))^2+(ps(3)-ps(1))^2));%计算写入mises应力到mpofset

end

2.5　节点应力插值计算

Abaqus插值的阶数与单元类型有关。一般来说，如果单元边中点处没有节点则为一阶插值，否则为二阶插值，相应的单元分别称为一阶或二阶单元，由于插值方法不同结果会略有不同。C3D8R单元为一阶单元，单元应力向节点应力的外插可采用线性插值，本文采用公用节点单元应力求和再平均的方法插值求得节点应力[1,3]。

算法3节点应力插值外推算法

load(′D:elmementofset.txt′) ; %载入几何集单元信息

load(′D: odesofset.txt′) ;%载入几何集节点信息

elmementofset(:,1)=[];

%置空单元编号列，避免单元编号与节点编号重复

nodestress=[];

%定义节点应力变量，存贮集合节点编号、坐标及节点应力

fori=1:1:length(nodesofset)%遍历几何集所属所有节点

[k,n]=find(elmementofset==nodesofset(i,1));

%在单元所属节点组列表中查找公用节点位置

sum=0;

sum2=0;

forij=1:1:length(k)%遍历共用节点

sum=sum+mpofset(k(ij,1),5);

%累加共用节点所在单元最大主应力

sum2=sum2+mpofset(k(ij,1),6);

%累加共用节点所在单元Mises应力

end

nodestress=[nodestress;nodesofset(i,1:4)sum/length(k)sum2/length(k)];

%平均求得节点应力，并与对应节点编号，坐标一并写入

End

3Abaqus节点应力内核调用输出法

Abaqus节点应力重构计算输出法体现了有限元软件节点应力输出的基本原理和过程，但对于一些复杂单元应力重构计算和插值外推算法十分复杂，实现较为繁琐。在Abaqus CAE的visualization模块，提供创建节点列表(nodelist)型路径(Path)，可以创建路径(Path)型XYData数据表，但节点列表的节点编号(label)，需要通过视图点选或手工输出，对于批量输出节点应力，工作量巨大，操作困难，且较易失误。此外，XYData数据表X坐标并不反映节点坐标，节点编号，数据解释性不好。若能通过Python编程生成批量输出的节点列表，则可利用Abaqus节点输出内核，避免Abaqus节点应力重构计算输出法单元应力重构和节点应力插值计算的麻烦。

3.1　Abaqus Session视图对象数据结构

Abaqus Session视图对象数据结构如图4所示，session对象主要包含图像、动画输出选项，视图定义，单元节点信息查询，数据表生成输出等功能。其主要功能都可在AbaqusCAE的Visualization模块实现，因此，利用Session视图对象进行Abaqus二次开发的研究较少。

图4　SESSION对象数据结构

3.2　Abaqus节点应力内核调用输出法

Abaqus节点应力内核调用输出法主要通过Python编程创建节点列表，利用session创建nodelist型path，其方法为session.Path()；创建path型xyDataObject批量输出节点应力，其方法为session.XYDataFromPath()。

算法4Abaqus节点应力内核调用输出算法

#本程序需打开Job-31.odb，在visualization模块相应输出界面运行。

from odbAccess import *

from abaqusConstants import*

myodb=openOdb(′Job-31.odb′)

nd=mdb.models[′Model-1′].parts[′Part-1′].sets[′Set-1′].nodes

v=1

nodelist=[]#定义节点列表nodelist为空列表

while v<=len(nd):

#遍历几何集Set-1所有节点，将节点编号写入节点列表nodlist

nodeslabel=mdb.models[′Model-1′].parts[′Part-1′].sets[′Set-1′].nodes[v-1].label

nodelist.append(eval(str(nodeslabel)))

v=v+1

else:

session.Path(name=′Path-1′,type=NODE_LIST,expression=((′PART-1-1′,(nodelist)),))

#创建以nodelist为节点序列的路径Path-1

pth = session.paths[′Path-1′]

session.XYDataFromPath(name=′XYData-1′,path=pth,includeIntersections=False,shape=DEFORMED,labelType=TRUE_DISTANCE)

#创建以Path-1为路径的XYData列表

print ′ok4′

xy1 = session.xyDataObjects[′XYData-1′]

#把XYData列表数据复制给变量xy1

with open(′D:/stress.txt′,′w′) as f4:

#把节点编号、坐标、应力值写入stress.txt

i=1

while i<=len(xy1):

index=nodelist[i-1]

stress=xy1[i-1]

nodescoors=mdb.models[′Model-1′].parts[′Part-1′].sets[′Set-1′].nodes[i-1].coordinates

f4.writelines(str(index)+′′)

f4.writelines(str(nodescoors)+ ′′)

f4.writelines(str(stress)+′ ′)

i=i+1

else:

f4.close()

myodb.close()

print ′ok5′

Stress.txt文件包含节点编号、坐标、应力值数据物理意义明确，但格式仍需通过Excel处理，方可提交Matlab进行分析处理。

4结果比较

Set-1共包含3284个节点，图5是Set-1所有节点应力着色面和等高线图，着色面图代表(X，Z)点对应应力大小，底面曲线为应力等高线，该图非常直观地反映了传感器膜片应力分布情况，且物理意义明确。

仅选如图6所示一组路径节点比较三种方法输出节点应力值。

图5　Set-1所有节点应力着色面和等高线

图6　一组随机选定的节点列表路径

结果表1所示，本文提出的Abaqus节点应力重构计算批量输出法和Abaqus节点应力内核调用批量输出法两种应力批量出方法，可以准确的定制输出节点应力，解决了Abaqus批量出节点应力的难题，且输出结果精度更高，意义明确，解释性好，为仿真结果的比较分析,特别是自动分析提供了便利。

表1　节点应力输出结果比较

5结论

基于Python-Matlab联合编程的Abaqus高级后处理技术可以便捷访问和处理计算结果，为计算结果后处理和自动分析提供了有效途径。基于重构计算和内核调用的节点输出算法均能有效实现节点应力的定制输出，且正确有效，精度较高，灵活性好，解释性好。重构计算输出法体现了有限元软件节点应力输出的基本原理和过程，但实现较为繁琐；内核调用输出法，仅从软件结构入手，不涉及节点应力输出的基本原理和过程，但简捷易行。

参考文献：

[1]ABAQUS Inc.ABAQUS Analysis User’s Manual/6.12.1 ABAQUS/Aqua analysis[K].2010.

[2]ABAQUS Inc.ABAQUS Scripting User’s Manual/6.12.1 ABAQUS/Aqua analysis[K].2010.

[3]曹金凤,石亦平.ABAQUS 有限元分析常见问题解答[M].北京：机械工业出版社,2009.

[4]曹金凤，王旭春，孔亮.Python语言在ABAQUS 中的应用[M]. 北京：机械工业出版社,2011.

[5]王勖成.有限单元法[M].北京：清华大学出版社，2003.

[6]张少实，庄茁.复合材料与粘弹力学[M].北京：机械工业出版社,2005.

[7]钟同圣，卫丰，王鸷. Python语言和ABAQUS前处理二次开发[J]. 郑州大学学报：理学版，2006，38(1):60- 64.

[8]张强，马永，李四超.基于Python的Abaqus二次开发方法与应用[J].舰船电子工程，2011，31(2):131-134.

[9]李猛,于存贵,崔二巍,等.ABAQUS 二次开发及在火炮参数化建模的应用[J].四川兵工学报,2013,34(9):41-43.

[10]徐凤军,高跃飞,常德顺,等.基于多体分析的火炮典型机构设计系统开发[J].四川兵工学报,2014,35(4):45-48.

(责任编辑唐定国)

收稿日期：2014-11-25

基金项目：国防预研基金(51328050101)

作者简介：任海峰(1978—)，男，博士研究生，讲师，主要从事火箭发动机健康监测技术研究；通讯作者：高鸣(1957—)，男，教授，博士，博士生导师，主要从事火箭发动机技术研究。

doi:10.11809/scbgxb2015.07.034

中图分类号：V435.1

文献标识码：A

文章编号：1006-0707(2015)07-0133-06

本文引用格式：任海峰，高鸣.Python-Matlab联合编程Abaqus高级后处理技术研究[J].四川兵工学报，2015(7):133-138.

Citationformat:RENHai-feng,GAOMing.ResearchonAbaqusAdvancedPostProcessingTechnologyBasedonPython-MatlabProgramming[J].JournalofSichuanOrdnance,2015(7):133-138.

ResearchonAbaqusAdvancedPostProcessingTechnology
BasedonPython-MatlabProgramming

RENHai-feng,GAOMing

(NavalAeronauticalandAstronauticalUniversity,Yantai264001,China)

Abstract:The further applications of finite element analysis put forward high request to flexibility and customization mode of result output. For it is difficult to output node stress in batches with Abaqus/CAE, the customization output methods were proposed utilizing programs of Python-Matlab by reconfiguration calculation and calling kernel of Abaqus. The mass customization output of node stress for the simulation results of bond stress sensor using in health monitoring of solid rocket motor were realized by proposed methods. The output results indicate that the suggested methods are correct, quite effective, higher precision, good flexibility and understandability.

Key words:solid rocket motor; bond stress sensor; Python-Matlab Programming; output in batches; advanced post processing of Abaqus; reconfiguration calculation; kernel call

【机械制造与检测技术】