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

【機械制造與檢測技術】