基于Abaqus/Python的三维复合材料层压板有限元建模参数化开发
伊鹏跃，于哲峰，汪海
（上海交通大学 航空航天学院，上海 200240）
摘要：讨论了Abaqus对象模型及网格对象的调用方法。利用Python脚本，开发了可应用连续壳或实体单元的三维复合材料层压板参数化建模程序。通过了人机交互，可选择层间是否使用界面单元，材料是否预制损伤。实例表明，该程序既可完成单个层压板及损伤的参数化三维建模，也可应用于装配形成的多板连接结构的建模，适用范围广泛。 关键词：Abaqus/Python；三维；复合材料；层压板；参数化建模 中图分类号：TP319;TB332
Parametric three-dimensional modeling of composite laminates based on
Abaqus/Python
Yi Pengyue, Yu Zhefeng, Wang Hai (School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai 200240) Abstract: Abaqus object model and invoking method of mesh object were discued.A parametric modeling program for three dimensional composite laminate was developed in which continuum shell or solid element is used.Human-computer interaction was realized by which laminate with interface and damage or not can be selected.Examples show that the program is able to be used not only in the modeling of one laminate and damage but also in the structures consisting of multiple plates by aembling.The application region of the program is wide.Keywords: Abaqus/Python; three dimension; composite; laminate; parametric modeling
引言
由于复合材料高比强度、比刚度及较强可设计性等特点在航空航天以及民用结构中的用量越来越大[1]。层压板式纤维增强复合材料在结构中应用广泛，其力学性能成为研究重点[2]。由于材料各向异性，层间性能差，损伤机理复杂，凡考虑层压板厚度、厚度方向应力分布或层间作用时，基本都需要使用三维模型。Abaqus具有丰富的材料本构和单元类型，并具备强大的二次开发能力，是复合材料结构分析的有力工具。
[3]基于Abaqus有限元分析软件，贾建东等通过不含界面层的三维复合材料层压板模型利用Umat编写材料本构，研究了层压板的冲击后剩余压缩强度；B.G.Falzon等[4]等基于三维
[5]实体模型研究了层压板的层内渐进损伤并进行了算例验证；顾亦磊等利用三维层压板模型
[6]使用Usdfld程序研究了层压板螺栓连接失效；Volnei Tita等利用三维模型对层压板静压痕试验进行了有限元仿真研究；张丽等[7]建立了含界面层的层压板三维有限元模型，利用Vumat子程序进行了复合材料层压板低速冲击作用下损伤分析；滕锦等[8]通过含粘性层的三维层压板模型研究了增韧复合材料的冲击损伤过程，朱炜垚等[9]仅在铺层角度不同的层间应用界面
[10-11]层简化三维模型，研究了层压板低速冲击下分层损伤；Craven等建立了含多处分层损伤与纤维损伤的层压板三维有限元模型，研究了带损伤层压板的拉伸、屈曲及后屈曲力学性能；王跃全等[12]利用层压板修补的三维有限元模型，进行了胶接贴补修理渐进损伤分析；张彦等[13]基于三维复合材料层压板模型对双悬臂梁试验进行仿真，验证了分层的临界能量释放率。
通过Abaqus的前处理模块进行三维层压板有限元建模是一个相当重复繁琐且耗时耗力的工作。如果开发参数化建模程序则可大大提高建模效率，促进层压板材料的性能分析和研究。本文采用Python脚本语言，基于孤立网格及单元节点信息，通过查找相应单元进行网格偏置，建立三维有限元网格，并赋予单元类型、材料属性与铺层角度等。开发完成了可包含界面层及可预制损伤的复合材料层压板三维参数化建模程序，应用了人机交互界面，对Abaqus对象模型、主要模块程序流程及实现方法进行介绍，并通过层压板三维建模例子演示程序的功能。
1 层压板建模与python脚本
三维复合材料层压板的建模繁琐，以基于网格偏置的建模方法为例，首先建立单层的几何模型，划分网格后将部件转化为孤立网格部件，进入mesh模块，通过Mesh edit功能对单层网格进行偏置，生成界面层或第二层复合材料单层。若为界面层，厚度值赋零，根据建模要求逐层偏置生成复合材料层压板三维模型，由于网格部件不具备几何信息，需为每一界面层与复合材料单层建立集合，利用集合为相应单元赋单元类型，材料属性，铺层角度等。如果铺层层数较多，通过用户界面进行操作，是一个重复工作量大且耗时耗力的过程。
Python脚本接口是Abaqus的二次开发接口之一，它扩充了Python的对象模型和数据类型。一般情况下，脚本接口主要用于前处理、后处理、自定义模块等。若前处理通过用户界面建模，需要大量的手动操作,而Python脚本使用区区几十行代码则可实现上百次操作的效果，若以此开发复合材料层压板参数化建模程序，可大大减少工作量，提高效率。
Abaqus/python主要有Seion对象、Mdb对象和Odb对象，本文主要对Mdb模型对象进行搜索调用。Mdb对象保存于模型数据库中，调用前需使用Import语句导入，图1为Mdb对象结构。通过对part对象的查询，给出了程序需用到的相关子对象，其中elements下的label、connectivity分别为单元编号、对应的节点编号，nodes下的coordinates、label分别为节点坐标、节点编号。
图1 Mdb对象结构 Fig.1 Mdb object model tree 2 含损伤的三维层压板参数化建模
本程序包括层压板的参数化建模部分和损伤建模部分，流程分别如图2(a)和图2(b)所示。 对于层压板的参数化建模，首先根据面内对网格的要求，例如是否需要局部加密等，可通过程序选择提前划分网格的平面部件，或者通过程序基于草图模块建立网格均匀的平面部件。网格偏置建模，需将平面部件转化为网格部件，然后根据层间是否引入界面层建立模型，包括不含界面层的层压板模型、各层间都有界面层的模型以及指定层间加入界面层的模型。由于模型中只含网格信息，为方便相应对象的选择，需要通过搜索单元编号为各层包括界面层建立单元集合。建立材料后，基于单元集合对相应的区域赋单元类型，对于铺层单元，可选择使用连续壳单元或实体单元，界面层单元使用Cohesive单元，对于不同的分析步类型，可分别选择Standard和Explicit单元，然后建立对应的截面属性并赋予相应的单元区域，最后根据提示对铺层角度进行赋值，至此三维层压板的建模完成。
对于损伤预制部分，建立损伤材料及截面属性后，会提示预制损伤的位置，可选择某一铺层或者界面层，然后需选择损伤区域形状，通过对该层所有单元的节点坐标进行遍历，并判断单元是否在损伤区域内，如此将该层所有单元划分为损伤单元集合与未损伤单元集合，为新建的损伤单元集合赋相应的截面属性。程序包含了两种基本的损伤形状，矩形与圆形，可通过循环对同一层单元集合多次预制损伤，取损伤形状的并集得到诸如椭圆形或花生壳状的复杂损伤形状。通过程序的循环，可建立多层、多处损伤。
至此，程序运行完成，用户只需要修改材料选项与参数，添加边界条件与载荷，即可提交计算。
(a)层压板参数化建模 (b)损伤预制
图2 程序流程图 Fig.2 Program flow chart 3 程序介绍
在脚本语言中首先利用from…import…语句导入相应模块，从而实现相关函数或对象的调用。程序运行时，为满足用户要求，需要通过交互界面进行选择或者输入，通过getInput()函数实现，如图3所示。
图3 交互界面 Fig.3 Interactive interface 描述层压板面内形状的平面部件类型属于3D的Shell部件，如图4所示，划分网格后，转化生成网格部件，然后输入铺层单层厚度，float()函数用于数据类型转换，通过len( )函数对网格部件单元数进行查询，选择所有shell单元的上表面，通过generateMeshByOffset( )命令生成实体单元即铺层首层，并激活deleteBaseElements选项删除壳单元，最后对单元重新编号。
##平面部件生成网格部件并偏置生成铺层首层 p = mdb.models['Model-1'].parts[planarpartname] p.PartFromMesh(name='Composite_laminates') input2=getInput('Please enter the PLY thickne(mm)','')#输入单层厚度 plyThickne=float(input2) p = mdb.models['Model-1'].parts['Composite_laminates'] e=p.elements plyelementNumber=len(e)#查询网格部件单元数
side2Elements = e[0:plyelementNumber]#选择要被偏置的shell单元 #网格的偏置操作
p.generateMeshByOffset(region=regionToolset.Region(
side2Elements=side2Elements), meshType=SOLID, totalThickne=plyThickne,
numLayers=1, offsetDirection=INWARD, shareNodes=True,
deleteBaseElements=True) p.renumberElement(startLabel=1, increment=1)# 单元重新编号
生成首层单元后需要继续偏置建立界面层和各铺层单元，以指定层间加入界面层的模型为例。首先通过循环结构输入需要引入界面层的层间编号，并记录到数组中。首先判断建立下一层铺层时是否需要建立界面层，若是，先偏置生成界面层单元并建立集合；若否，则直接生成下一层铺层单元并建立集合；如此循环直到偏置完成。由于各单层单元总数相同，单元的编号按照生成的先后顺序从小到大排序，所以偏置操作的基础单元可通过计算上次生成的单元编号，根据对象结构进行选择，单元集合的建立同样通过计算相应单元的编号进行选择。
##指定层间加入界面层的三维层压板网格偏执建模
intface_insert=1
intface_location=[] #循环输入需引入界面层的层间编号
while intface_insert:
input7=getInput('Enter interface location (e.g Interface3:3)','')
intf_loc=int(input7)
intface_location.append(intf_loc)
input8=getInput('More interface or not (Yes:1,No:0)','')
intface_insert=int(input8)
intface_locationLength=len(intface_location)
intface_locationIndex=0
for x in range(1,plyNumber,1):#循环逐层建模
if x==intface_location[intface_locationIndex] and intface_locationIndex
#网格偏置生成界面层
p = mdb.models['Model-1'].parts['Composite_laminates']
f = p.elements
#计算单元编号选择偏置操作的基础单元
face2Elements = f[plyelementNumber*(intface_locationIndex+x-1):(intface_locationIndex+x)*plyelementNumber]
#网格偏置生成下一层
p.generateMeshByOffset(…)
#界面层单元集合创建
p = mdb.models['Model-1'].parts['Composite_laminates']
e = p.elements
#计算单元编号选择相应单元
elements = e[(intface_locationIndex+x)*plyelementNumber:(intface_locationIndex+x+1)*plyelementNumber]
p.Set(elements=elements, name='Interface%d'%(x))#为相应单元建立集合
intface_locationIndex=intface_locationIndex+1
#网格偏置生成铺层
p = mdb.models['Model-1'].parts['Composite_laminates']
f = p.elements
face2Elements = f[(intface_locationIndex+x-1)*plyelementNumber:(intface_locationIndex+x)*plyelementNumber]
p.generateMeshByOffset(…)
#铺层单元集合建立
p = mdb.models['Model-1'].parts['Composite_laminates']
e = p.elements
elements = e[(intface_locationIndex+x)*plyelementNumber:(intface_locationIndex+x+1)*plyelementNumber]
p.Set(elements=elements, name='PLY%d'%(x+1))
网格偏置建模完成后，选择单元类型，建立材料及相应截面属性，通过循环基于单元集合为对应区域赋属性，如图5所示。同理对铺层角度赋值，选择之前建立的铺层单元集合，利用orientation功能逐层定义材料角度。
##铺层角度赋值
for x in range(1,plyNumber+1,1):#对各铺层循环
input3=getInput('Please enter the angle of the PLY%d'%x,'')#输入当前铺层的角度
plyAngle=int(input3)
p = mdb.models['Model-1'].parts['Composite_laminates']
region = p.sets['PLY%d'%(x)]#选择相应铺层单元集合
orientation = mdb.models['Model-1'].parts['Composite_laminates'].datums[coordinate_id]#选择参考坐标系 #为相应铺层赋角度
mdb.models['Model-1'].parts['Composite_laminates'].MaterialOrientation(region=region,
orientationType=SYSTEM, axis=AXIS_3, localCsys=orientation, fieldName='',
additionalRotationType=ROTATION_ANGLE, additionalRotationField='',
angle=plyAngle, stackDirection=STACK_3)
对于损伤区域的建立，首先遍历各未损伤单元集合，然后根据对象模型树中数据结构对集合中的单元信息进行查询判断。以圆形损伤区域为例，对集合中的每一单元的节点进行遍历，若该单元所有节点处于区域内，将其编号添加到损伤单元对应的数组中，如果某一节点不在该区域内，则将其编号添加到另一数组中，如此划分开来，遍历结束后根据两组数组中的单元编号分别建立损伤单元集合与未损伤的单元集合。程序运行结束后建立的各铺层和界面层及损伤区域的单元集合，如图6所示，根据损伤区域的单元类型建立对应的损伤截面属性，对损伤的铺层材料重新赋角度。
##损伤单元集合的建立
for element in eSET:#对该集合中单元进行遍历
state=0
nE=len(element.connectivity)
#对单元对应的节点进行遍历并通过节点坐标判断单元是否位于损伤区域内
for n in range(0,nE):#节点编号遍历
x=allnodes[element.connectivity[n]].coordinates[0]#当前节点x轴坐标查询
y=allnodes[element.connectivity[n]].coordinates[1] #当前节点y轴坐标查询
arcdistance=sqrt((x-coord_xc)**2+(y-coord_yc)**2)
if arcdistance>r+tolerance:#当前节点是否在圆形损伤区域内
state=1
j.append(element.label)#将当前节点对应的单元的编号添加到未损伤集合中
break
if state==0:
i.append(element.label) #将当前节点对应的单元的编号添加到损伤集合中
图4平面部件
图5 截面属性与铺层角度 图6 单元集合 Fig.4 Planar part
Fig.5 Section and ply angle
Fig.6 Element sets
4 应用实例
在圆形复合材料层压板冲击模型中，需要对冲击点处网格进行细化，基于预先划分网格的平面部件，运行程序得到的层压板模型如图7所示。在研究材料损伤对受压层压板应变场的影响时，需建立花生壳状的分层区域，运行程序，完成层压板建模后，选择预制损伤，损伤区域通过取圆心不同的两个圆形区域的并集得到，如图8所示。
图7 三维复合材料层压板模型 图8 花生壳状分层
Fig.7 Three dimensional composite laminates model
Fig.8 Peanut shape delamination
在复合材料材料层压板螺栓连接孔边挤压破坏的研究中，基于预先划分网格的平面部件，运行程序得到层压板模型，并将其与金属钉模型进行装配得到钉连接模型，如图9(a)所示。在复合材料层压板的双搭接修补模型拉伸强度分析中，运行程序得到母片与补片的层压板模 形，并在母片的修补区域通过偏置生成胶层模型，然后进行双搭接的装配并利用Merge命令合并重节点从而将结构转化为整体，模型如图9(b)所示。
(a)
图9 复合材料层压板螺栓连接与双搭接
Fig.9 Bolted Joint and double lap joint of composite laminates
(b)
5 结论
(1)构建程序流程，嵌入内部循环，加入交互界面，既保证了程序结构的简明易实现，又满足了参数输入和循环操作的要求。
(2)分析了对象模型的数据结构，并应用到单元与节点的搜索筛选，实现了基于网格信息的参数化建模过程，为基于网格信息的参数化操作提供了参考。 (3)根据有限元分析要求，程序可建立多种形式的三维层压板模型,并可实现损伤建模，也可应用于基于层压板的多板结构的建模，适用范围广，建模效率高。 参考文献：
(a)
[1]张小娟,张博平,张金奎 等.基于凹坑深度的复合材料低速冲击损伤分析[J].实验力学, 2010,25(3):234-238.
[2]矫桂琼，贾普荣.复合材料力学[M].西安：西北工业大学出版社，2008.3.[3]贾建东，丁运亮，刘晓明.复合材料层合板冲击后剩余强度的工程估算方法和有限元模拟分析[J].南京航空航天大学学报，2010，42(3),335-339.[4]Falzon B G, Apruzzese P.Numerical analysis of intralaminar failure mechanisms in composite structure[J].Composite Structures,2011,93(2):1039-1053.[5]顾亦磊，赵美英.复合材料层合板螺栓连接失效分析[J]航空计算技术，2006,36,(2):110-113.[6]Tita V, Carvalho J, Vandepitte D.Failure analysis of low velocity impact on thin composites laminates [J].Composites Structure,2008,83:413-428.[7]张丽,李亚智,张金奎.复合材料层合板在低速冲击作用下的损伤分析[J].科学技术与工程,2010,5(10):1170-1174.[8]滕锦，李斌太，庄茁.z-pin增韧复合材料层合板低速冲击损伤过程研究[J].工程力学，2006,23(1)：209-216.[9]朱炜垚,许希武.复合材料层合板低速冲击损伤的有限元模拟[J].复合材料学报,2010,6(27):200-207.[10]Craven R, Sztefek P, Olon R.Investigation of impact damage in multi-directional tape laminates and its effect on local tensile stiffne[J].Composites Science and Technology,2008,68:2518–2525.[11]Craven R, Iannucci L, Olon R.Delamination buckling: A finite element study with realistic delamination shapes, multiple delaminations and fibre fracture cracks[J].Composites: Part A,2010,41:684-692.[12]王跃全,童明波,朱书华.复合材料层合板胶接贴补修理渐进损伤分析[J].复合材料学报,2011,28(3):197-202.[13]张彦.纤维增强复合材料层合结构冲击损伤预测研究[D].上海: 上海交通大学博士学位论文,2007.
作者简介：伊鹏跃(1988-),男，山东泰安人，硕士。研究方向：飞行器结构强度。E-mail：yipengyue@sjtu.edu.cn
