利用Python提取ABAQUS的计算结果（ODB）信息

一个金属长方体，我们需要对其做拉伸的加载约束示意图如图1，并在完成后采用Python命令流读取参考点的位移、体积、应变随加载时间的变化情况。

 

                                                 图1 金属长方体约束加载示意图 

 

01

Python库导入声明

 

对于图上的命令流，我们分为几部分拆开研究，首先是下面第一部分需要在Abaqus导入相关Python库的代码，比如matplotlib库、numpy库、xlwt库和math库。

 

#! /user/bin/python
# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from visualization import *
from odbAccess import *
import xlwt
import math

 

• 第3行导入Python经常用于绘图的matplotlib库；

• 第4行导入Python处理数组矩阵的numpy库；

• 第5行和第6行用于导入Abaqus的后处理功能，即visualization和odbAccess；

• 第7行用于导入Python与excel进行交互的xlwt库；

• 第8行用于导入Python的math库。

在利用Python生成excel数据之前，先采用一个getInputs函数（代码如下）生成图2所示的对话框与用户进行交互，需要用户在在其中输入模型的名字（model name），部件例子的名字（instance name）以及后处理odb的名字（odbname）。这样做有个好处，每次都可以根据不同模型、不同部件实例和相应的后处理odb名字进行相应结果提取。

 

modelName,instanceName,odbname= getInputs( 
 fields=(('Model Name:', 'Test'),('Instance Name:', 'Part-2-1'),('odbname:', 'TEST.odb')), 
  label='Enter information',  dialogTitle='Enter information.')

 

                                                    图2 getInputs函数生成的对话框

 

02

计算初始边长

 

 

在计算体积和各边应变之前，需要计算长方体不同边长即长、宽、高的初始长度。代码如下：

node = mdb.models[modelName].rootAssembly.instances[instanceName].nodes
Xmin = 9999
Xmax = -9999
Ymin = 9999
Ymax = -9999
Zmin = 9999
Zmax = -9999
for i in range(len(node)):   
    x = node[i].coordinates[0]    
    y = node[i].coordinates[1]    
    z = node[i].coordinates[2]    
    if Xmin > x:        
        Xmin = x    
    elif Xmax < x:        
        Xmax = x            
    if Ymin > y:        
        Ymin = y    
    elif Ymax < y:        
        Ymax = y            
    if Zmin > z:        
        Zmin = z        
        continue    
    if Zmax < z:        
        Zmax = z        
        continue
print 'Xmin,Xmax,Ymin,Ymax,Zmin,Zmax=',(Xmin, Xmax, Ymin, Ymax, Zmin, Zmax)
CubeWidth=abs(Xmax-Xmin)  
# the length of model along x directionCubeHeight=abs(Ymax-Ymin) 
# the length of model along y directionCubeLength=abs(Zmax-Zmin) 
# the length of model along z direction

• 第1行基于之前输入的modelname，instancename来调用部件实例（instance）的节点。

• 第8行for i in range(len(node))：表示对每个节点进行循环，每次循环中都要提取该节点的x，y，z坐标信息。

• Xmin，Xmax为通过循环找出X方向的最小、最大X坐标，Ymin、Ymax以及Zmin、Zmax都同理。计算这部分的代码可能大家表面无法直接看懂理解，我举一个例子好了。

 

例子：对于第12行到第15行，比如当i=0时，随便提取了第1个节点的x，y，z信息（0,0.1,0.2），由于Xmin=9999>x=0，所以x=0赋予Xmin；同样的，当i=1时，又随便提取模型的一个节点(0.05,0.1,0.2)，由于现在Xmin=0<x=0.05，则执行Xmax=-9999<x=0.05，所以x=0.05赋予Xmax，而Xmin不变；当i=2时，随便提取一个节点(0.07,0.1,0.2)，注意到此时，Xmin=0<x=0.07并且Xmax=0.05<x=0.07，因此x=0.07赋予Xmax，即此时的最大值赋予给Xmax，而Xmin不变；而如果i循环到节点 (0.1,0.1,0.2)，注意到此时，Xmin=0<x=0.1并且Xmax=0.05<x=0.1，因此x=0.1赋予Xmax，即此时的最大值赋予给Xmax，而Xmin不变；由于x=0是模型x方向的最小值，x=0.1是模型x方向的最大值，因此如果当i循环到又一个节点(0.06,0.1,0.2)，此时，Xmin=0<x=0.06并且Xmax=0.1>x=0.06，则Xmin和Xmax则处于不变，因为x在这两者之间，因此通过以上程序就可以找出模型在x方向的最大值和最小值（结果见图3），然后最大值减去最小值则得到了模型沿x方向的长度。同样的方式得到模型沿y方向和z方向的长度。

 

 

图3 查看Xmin、Xmax、Ymin、Ymax、Zmin、Zmax

 

另一方面，我们可以如图4进行相应验证，我们先查看下模型在x方向的宽度，在Abaqus中的query中选取distance，然后选中模型中的两个点，可以看到这两个点的坐标是point1（0,0.1,0.2），point2（0.1,0.1,0.2），则两点之间在x方向的距离为0.1，即为模型的宽度（模型沿x方向的长度）。Abaqus模型计算的结果与之前程序计算的结果相互验证了。

 

图4 查看模型在x方向的宽度

 

 

03

 

 

用python将位移导入excel

odb=openOdb(path=odbname)
wbkName='NodalDisplacement'
wbk1=xlwt.Workbook()
sheet=wbk1.add_sheet('Sheet1')
RefPointSet=odb.rootAssembly.nodeSets['U-XYZ']
frameRepository=odb.steps['Step-1'].frames
for i in range(3):    
    sheet.write(0,i,('U'+str(i+1)+'(U-XYZ)'))

• 第1行表示打开odb计算结果文件。

• 第2行表示生成的excel文件名为NodalDisplacement。

• 第3行和第4行表示添加Sheet1工作表。

• 第5行建立节点SetU-XYZ，方便后续对该Set进行位移提取。

• 第6行表示Step-1里的帧(frames)对象。

• 第7行和第8行通过循环在excel第1行第i列（也就是代码中的0，i），写入U1(U-XYZ)，U2(U-XYZ)，U3(U-XYZ)，如图5所示，其中Python编程用+可以实现字符连缀，str(i)可以将数字变量i改为字符i。

iframes=len(frameRepository)
RefDataU1= np.zeros((iframes,1))
RefDataU2= np.zeros((iframes,1))
RefDataU3= np.zeros((iframes,1))
jingdu=15
for i in range(iframes):    
    U=frameRepository[i].fieldOutputs['U']    
    RefU=U.getSubset(region=RefPointSet)   
    RefUValues=RefU.values    
    RefDataU1[i]=RefUValues[0].data[0]   
    sheet.write(i+1,0,round(RefDataU1[i],jingdu))   
    RefDataU2[i]=RefUValues[0].data[1]    
    sheet.write(i+1,1,round(RefDataU2[i],jingdu)) 
    RefDataU3[i]=RefUValues[0].data[2]   
    sheet.write(i+1,2,round(RefDataU3[i],jingdu))
wbk1.save(wbkName+'.xls')

• 第一行是统计frames的帧数。

• 第2行到第4行基于numpy定义三个数组存储三个方向的位移，并事先赋予初值零。

• 第5行设置一个变量，方便用于控制输出到excel里数的精度。

• 第6行到第15行就是通过for循环，一帧一帧地输出U-XYZ这个点Set的三个方向的位移值，并存储到excel。U-XYZ这个点Set如图5所示。

图5 U-XYZ点Set

 

• 对于第7行表示提取该帧的位移U的场变量并赋值给U，然后在第8行采用getSubset基于U提取RefPointSet区域的位移值。

• 图6在ABAQUS命令行接口（基于Python的Abaqus静力分析操作实例）表示的应该是最后一次循环后各变量的值，可以看出RefU表示的nodeSets['U-XYZ']的位移值，图6中显示为一个索引；而RefUValues则表示该节点集中每个节点的位移值，体现在values中；而values[0]表示第一个节点的位移值，由于存在三个方向的位移，U1,U2,U3,则后续通过data[0]，data[1]，以及data[2]分别表示三个方向的位移值。

 

图6 运行相关结果查看（通过ABAQUS命令行接口）

 

特别的，由于U-XYZ只有一个点存在于Set中，如果我们在这里用values[1]，看会出现什么样的效果？如图7。可以明显看到Sequence index out of range表示超出了范围，因为RefUValues[1]表示提取第二个节点，但我们这里只有一个节点存在于U-XYZ点Set中，所以出现了范围超过的报错提示。

 

图7 可以存在RefUValues[1]吗？（通过ABAQUS命令行接口）

 

需要说一下的，如果大家觉得ABAQUS命令行接口的查看范围太窄，可以往上拉动窗口，这样就可以看到更多的代码，如图8所示，但与之相反的就是模型就会显得很小，与我们一般的常规abaqus操作有点相反。

图8 ABAQUS命令行接口显示更多的内容

 

再对代码中的Python round函数做个说明，通过下面的例子我们可以看到round函数是四舍五入函数，一般为四舍五入后为整数，round(a,b)中的a为需要进行四舍五入的数，而b为需要保留的小数位数，如图9所示。

 

图9 关于Python round 函数的例子

 

最后就是第17行代码，这个代码也很重要，不输入这个代码就不会生成excel表，一般生成excel位置为设置好的工作目录（比如D:\temp），生成好的excel如图10所示。而excel总共有625行，说明一共有625帧，也就是上述代码的变量iframes为625。

 

图10 生成的NodalDisplacement.xls的excel文件

04

结果验证

将上述excel表格中的结果绘制成曲线，即U-XYZ点三个方向的位移随帧数变化的曲线，如图11所示。根据图12，再结合odb文件中的变形图（第623帧）可以判断U-XYZ点在后期的X(U1)为负，Y(U2)为正，Z(U3)为正，与图11曲线是一致的。

 

图11 U-XYZ点三个方向的位移随帧数变化的曲线

 

图12 模型后处理结果图（第623帧）

 

最后再与在Abaqus里进行GUI（即点点点）对U-XYZ三个方向位移处理的结果图进行对比，即通过XY Data对U-XYZ三个方向的位移进行输出绘图，注意到坐标轴为时间，与Python后处理图11的帧数是一一对应的。可以知道，Python后处理结果与ABAQUS的GUI操作得到的结果几乎完全一致。

 

图13 Abaqus里进行GUI操作的后处理结果（关于U-XYZ三方向位移）

该模型关于体积、应变、应力等方面（包括场输出和历史输出）的Python导入excel的后处理留在下期，敬请关注。

最后，欢迎大家通过微信公众号联系我们。

微信公众号：320科技工作室。

同时，更多ABAQUS二次开发的精彩内容，欢迎大家关注微信公众号： 土木科研编程与数值模拟