基于ABAQUS平台的机电虚拟裂纹闭合法

周立明+孟广伟+李鹏+李锋+李宵琳

摘要：为提高求解含裂纹压电体的能量释放率的计算效率和精度以及减少程序编写与调试的工作量，针对四节点平面压电单元提出了机电虚拟裂纹闭合法计算公式。基于通用商业有限元软件ABAQUS，开发了哑节点断裂压电单元，编写用户自定义子程序UEL，该单元可独立求解总能量释放率分量，对不同材料和裂纹长度的含裂纹压电体的能量释放率分量进行了求解，与理论解作了对比，讨论了不同形式网格对求解精度的影响。数值算例结果表明，哑节点断裂压电单元具有精度高、简单方便和对网格尺寸不敏感等优点。

关键词：裂纹；机电虚拟裂纹闭合法；能量释放率；哑节点断裂压电单元； ABAQUS

中图分类号：TB115 文献标识码：A

越来越多的压电器件在苛刻环境条件下服役，裂纹是导致构件失效的主要因素，对构件进行断裂分析的第一步便是断裂参数的求解。目前计算断裂参数的方法＼[1-2＼]有外推法、J积分、扩展有限元法、虚拟裂纹扩展法和虚拟裂纹闭合法。外推法要求裂尖处具有特别细的网格或采用奇异元；J积分表达式繁琐，不容易被工程师广泛采用；扩展有限元法裂尖处单元需采用裂尖渐近位移场函数进行加强；虚拟裂纹扩展法需要两次有限元分析；虚拟裂纹闭合法对网格尺寸不敏感，表达式简单，容易编程，只需一步有限元分析，计算精度高，在解决实际工程问题中发挥着重要作用。

Rybicki和Kanninen＼[3＼]于1977年提出虚拟裂纹闭合技术，将该方法应用于求解含裂纹结构的应变能释放率。Raju＼[4＼]和Xie等＼[5＼]对虚拟裂纹闭合法进行了数学解释。Xie等＼[6-8＼]提出了哑节点断裂单元，为虚拟裂纹闭合法的发展作出了贡献。虚拟裂纹闭合法已应用于复合材料、功能梯度材料、加强结构的断裂分析与评估、岩土材料等领域＼[9-10＼]，但机电耦合场下求解结构断裂参数的虚拟裂纹闭合法还未见报道。

工程结构由于几何形状、材料属性和加载方式的复杂性，不得不依赖于数值方法。随着科学技术的发展，用数值方法求解断裂参数变得切实可行，很多计算方法被用来求解断裂参数，例如有限差分法＼[11＼]、无网格方法＼[12＼]、杂交元＼[13＼]、边界元＼[14＼]和光滑有限元＼[15＼]等，但由于缺少商业软件的支持，这些数值方法的工程实际应用相对缺乏。选用有限元软件ABAQUS为平台，可直接从软件计算的结果中提取相关信息，通过编写用户自定义单元（UEL）子程序，实现含裂纹压电材料断裂参数的计算，可极大地减少程序编写和调试的工作量，程序一旦得到验证，很方便应用到工程实际的断裂分析中，提高程序的通用性和计算效率。

湖南大学学报（自然科学版）2015年

第10期周立明等：基于ABAQUS平台的机电虚拟裂纹闭合法

在含裂纹压电材料中，裂纹扩展单位长度所需要的能量称为总能量释放率，包含机械能释放率GM和电能释放率GD，GM为裂纹尖端前方的应力在裂纹改变位移上所做的功，GD为裂纹尖端前方的电位移在裂纹变化后裂纹面电势差上所做的功。本文针对四节点平面压电单元提出了机电虚拟裂纹闭合法计算公式。以通用商业有限元软件ABAQUS为平台，开发了哑节点断裂压电单元，编写用户自定义子程序UEL，该单元可独立求解总能量释放率分量，对不同材料和裂纹长度的含裂纹压电体的总能量释放率分量进行了求解，并与理论解作了对比，讨论了不同形式网格对求解精度的影响。

1机电虚拟裂纹闭合法

针对求解含裂纹压电材料的总能量释放率的需要，提出了针对四节点平面压电单元的机电虚拟裂纹闭合法计算公式。考虑一含裂纹压电体，如图1所示，裂纹长度为a，厚度为B，总能释放率G为产生面积为ΔA的新裂纹面所需要的能量，Δa为裂纹扩展量，于是

G=GI+GII+GD；（1）

GI=limΔa→012BΔa∫Δa0σ1yyΔa-r，0Δv2r，πdr；（2）

GII=limΔa→012BΔa∫Δa0τ1yyΔa-r，0Δu2r，πdr；（3）

GD=limΔa→012BΔa∫Δa0D1yΔa-r，0Δφ2r，πdr。（4）

式中：σ1yy，τ1yy和D1y分别为原始裂纹尖端处的法向应力、切向应力和法向电位移；Δv（2），Δu（2）和Δφ（2）分别为裂纹虚拟扩展到a+Δa时裂纹面上的张开位移、相对滑动位移和电势差；GI和GII为机械能释放率分量；GD为电能释放率。

图1机电虚拟裂纹闭合法示意图

Fig。1Electromechanical virtual crack closure technique

基于势能的改变与将裂纹闭合一个扩展增量所需的功等效，提出机电虚拟裂纹闭合法。将虚拟裂纹闭合算法在其提出的假设基础上进行了横向扩展，即在虚拟裂纹扩展过程中计入电势以及位移的作用，相对应的压电单元，将电势作为一个“位移”分量进行考虑，式（2）-式（4）可改写为：

GI=yΔ2BΔa；（5）

GII=xΔ2BΔa；（6）

GD=Δ2BΔa。（7）

式中：x和y，Δ和Δ，，Δ分别为局部坐标系（，）下节点力分量、张开位移分量、节点电荷和电势差。

由式（5）-式（7）可知：1）机电虚拟裂纹闭合法可分别计算GI，GII，GD；2）能量释放率的计算仅仅包含节点力与节点位移、节点电荷与节点电势差，这些变量可从有限元软件中输出；3）避免了应力和电位移的积分，容易与有限元分析相结合。

2 哑节点断裂压电单元

基于ABAQUS平台，利用哑节点断裂压电单元来实现二维线状含裂纹压电断裂力学问题，通过编写用户自定义单元子程序UEL来实现。哑节点断裂压电单元的定义及其节点编号如图2所示，含5个节点。节点1和节点2对应于裂纹尖端，在其节点间放置有特殊刚度的弹簧，节点3和节点4位于裂纹尖端的后面，节点5在裂纹尖端的前面。在ABAQUS中，单元所具有的完整的节点矢量为：

U={u1，v1，1，u2，v2，2，u3，v3，3，u4，v4，

4，u5，v5，5}T。（8）

裂纹尖端处的节点力和节点电荷为：

Fx=Kx（u1

Symbolm@@ u2）， Fy=Ky（v1

Symbolm@@ v2），

Q=K（1

Symbolm@@ 2）。（9）

式中：u1和v1，u2和v2分别为节点1和节点2在整体坐标系（X，Y）下位移分量；1和2分别为节点1和节点2在整体坐标系（X，Y）下电势；Kx和Ky分别为力场下X和Y方向的弹簧刚度；K为电场的弹簧刚度。

节点3、节点4和节点5被用来从ABAQUS结果中提取相关信息，对单元的刚度矩阵并没有实际贡献，称为“哑节点”，该单元为哑节点断裂压电单元。虚拟裂纹扩展量是主节点1和哑节点5之间的距离：

Δa=x5-x12+y5-y12。（10）

式中：（x1，y1）和（x5，y5）分别为节点1和节点5在整体坐标系（X，Y）下的坐标。如果

SymbolDA@ a在每个增量步中都通过位移来更新，裂纹方向随之更新，在大变形分析问题时很有用。

图2哑节点断裂压电单元

Fig。2Fracture of piezoelectric element

with dummy nodes

应变能释放率必须在裂纹尖端处的局部坐标系（，）下计算，X轴和轴之间的夹角为：

cosθ=x5-x1Δa，sinθ=y5-y1Δa。（11）

通过简单的矢量投影关系，在局部坐标系（，）下节点力和节点电荷为：

x=Fxcosθ+Fysinθ；

y=-Fxsinθ+Fycosθ；

=Q。（12）

张开位移和电势差为：

Δ=Δucosθ+Δvsinθ；

Δ=-Δusinθ+Δvcosθ；

=Δφ。（13）

能量释放率可由式（12）和式（13）代入式（5）-式（7）计算得到。

3数值算例

如图3所示，一含中心裂纹的压电体，裂纹长度为2a，压电材料的极化方向为P，边长2l=60 cm，受均匀拉伸

SymbolsA@

SymboleB@ =1×105Pa和电位移D

SymboleB@ =7。4×10-5C/m2的作用，采用PZT4，P7和PZT5H 3种材料进行数值模拟，材料属性见表1，对于该裂纹问题，能量释放率的理论解为＼[16＼]。

PZT4：

G=（0。362 9

SymbolsA@

SymboleB@ 2+0。373

SymbolsA@

SymboleB@ E

SymboleB@

Symbolm@@ 138。3E

SymboleB@ 2）×

10-10a。（14）

P7：

G=（0。406 8

SymbolsA@

SymboleB@ 2

Symbolm@@ 0。446

SymbolsA@

SymboleB@ E

SymboleB@

Symbolm@@ 428。5E

SymboleB@ 2）×

10-10a。 （15）

PZT5H：

G=（0。424 8

SymbolsA@

SymboleB@ 2

Symbolm@@ 0。695 2

SymbolsA@

SymboleB@ E

SymboleB@

Symbolm@@ 389。44E

SymboleB@ 2）×

10-10a。（16）

其中E

SymboleB@ 与D

SymboleB@ 关系为：

D

SymboleB@ =c11c33-c13e31c11c33-c213σ

SymboleB@ +

d33+c33e231-2c13e31e33+c11e233c11c33-c213E

SymboleB@ 。（17）

图3含中心裂纹压电体模型

Fig。3Piezoelectric model with a centre crack

由于结构和载荷的对称性，取1/4结构进行求

解，顶部施加相应的应力和电位移，约束左端所有节点的水平方向位移和底部裂尖以右的节点电势（即底部电势为零）和竖直方向位移，为验证机电虚拟裂纹闭合法（EMVCCT）的可靠性，当裂纹长度2a=2 cm时，将结构离散为I型（单元：15×15），Ⅱ型（单元：30×30）和Ⅲ型（单元：60×60） 3种均匀分布网格形式，该结构为I型裂纹，仅需考虑含裂纹压电体的GI和GD，表2为EMVCCT计算得到的GI，GD和理论解。

由表2可以看出，EMVCCT在3种网格离散形式、3种材料下均得到了精度较高的GI和GD，与理论解误差最大不超过3。04%，3种网格所得精度基本一致，可见，该方法不仅具有较高的精度，且对网格的尺寸大小不敏感。

为方便计算不同裂纹长度下结构的GI和GD，采用Ⅲ型单元离散形式，表3给出了3种材料在不同裂纹长度下得到的GI和GD，并与理论解做了比较，从结果可以看出EMVCCT的计算结果比解析解得到的结果要小，是由于插值函数使用了“协调和完整的位移函数”，连续体离散后刚度会有所增加，求解值相对实际值要小，因此，EMVCCT的计算结果比解析解得到的结果要小；从结果还可以看出EMVCCT具有较高的精度，进一步验证了EMVCCT具有表达式简单，容易编程，只需一步有限元分析，计算精度高的优点。

4结论

本文针对四节点平面压电单元提出了机电虚拟裂纹闭合法，基于通用商业有限元软件ABAQUS平台，开发了哑节点断裂压电单元，编写用户自定义子程序UEL，以PZT4和P7、PZT5H的压电平板的中心裂纹问题为例，求解了不同网格离散形式和裂纹长度下结构的总能量释放率，并与理论解做了对比，结论如下：

1）该方法表达式简单，容易编程，只需一步有限元分析，对网格的尺寸大小不敏感，具有较高的计算精度。

2）该方法基于有限元软件ABAQUS可直接从软件计算的结果中提取相关信息，极大地减少了程序编写和调试的工作量。

参考文献

[1]LESKI A。 Implementation of the virtual crack closure technique in engineering FE calculations[J]。Finite Element in Analysis and Design，2007，43（3）：261-268。

[2]FRIES T P， BELYTSCHKO T。 The extended/ generalized finite element method： An overview of the method and its applications＼[J＼]。International Journal for Numerical Methods in Engineering， 2010， 84（3）：253-304。

[3]RYBICKI E F， KANNINEN M F。 A finite element calculation of stress intensity factors by a modified crack closure integral＼[J＼]。 Engineering Fracture Mechanics， 1977， 9： 931-938。

[4]RAJU I S。 Calculation of strainenergy release rates with highorder and singular finiteelements＼[J＼]。 Engineering Fracture Mechanics， 1987，28： 251-274。

[5]XIE D， BIGGERSJR S B。 Calculation of transient strain energy release rates under impact loading based on the virtual crack closure technique＼[J＼]。 International Journal of Impact Engineering， 2007， 34（6）： 1047-1060。

[6]XIE D， BIGGERSJR S B。 Progressive crack growth using interface element based on the virtual crack closure technique＼[J＼]。 Finite Elements in Analysis and Design， 2006， 42（11）：977-984。

[7]HE W， LIU J， XIE D。 Numerical study on fatigue crack growth at a webstiffener of ship structural details by an objectedoriented approach in conjunction with ABAQUS＼[J＼]。 Marine Structures， 2014， 35： 45-69。

[8]XIE D， SHERRILL B， BIGGERS J。 Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique。 Part II： Sensitivity study on modeling details＼[J＼]。 Engineering Fracture Mechanics， 2006， 73（6）： 786-801。

[9]SENTHIL K， AROCKIARAJAN A， PALANINATHAN R。 Defects in composite structures： Its effects and prediction methods-a comprehensive review ＼[J＼]。 Composite Structures，2013， 106： 139-149。

[10]周立明，孟广伟，王晖，等。 基于光滑有限元的含裂纹复合材料的虚拟裂纹闭合法＼[J＼]。 湖南大学学报：自然科学版，2014，41（8）：13-17。

ZHOU Liming， MENG Guangwei， WANG Hui，et al。 Virtual crack closure technique based on smoothed finite method for composite meterials with cracks＼[J＼]。 Journal of Hunan University：Natural Sciences，2014， 41（8）： 13-17。（In Chinese）

[11]LIAO D M， ZHANG B， ZHOU J X。 Using finite difference method to simulate casting thermal stress＼[J＼]。 China Foundry， 2011， 8（2）： 177-181。

[12]龙述尧， 张国虎。 基于MLPG法的动态断裂力学问题＼[J＼]。 湖南大学学报：自然科学版，2012， 39（11）：41-45。

LONG Shuyao， ZHANG Guohu。 An analysis of the dynamic fracture problem by the meshless local PetrovGalerkin method＼[J＼]。 Journal of Hunan University：Natural Sciences， 2012， 39（11）：41-45。 （In Chinese）

[13]平学成， 陈梦成， 谢基龙，等。 基于新型裂尖杂交元的压电材料断裂力学＼[J＼]。 力学学报， 2006， 38（3）：407-413。

PING Xuecheng， CHEN Mengcheng， XIE Jilong， et al。Fracture mechanics researches on piezoelectric materials based on a novel cracktip hybrid finite element method＼[J＼]。 Acta Mechanica Sinica， 2006， 38（3）：407-413。 （In Chinese）

[14]BENEDETTI I， ALIABADI M H， MILAZZO A。 A fast BEM for the analysis of damaged structures with bonded piezoelectric sensors＼[J＼]。 Computer Methods in Applied Mechanics and Engineering， 2010， 199（9）： 490-501。

[15]ZHOU L M， MENG G W， LI F， et al。 Cellbased smoothed finite element methodvirtual crack closure technique for a piezoelectric material of crack[J]。Mathematical Problems in Engineering，2015，371083：1-10。

[16]ZHANG T Y， QIAN C F， TONG P。 Linear electroelastic analysis of a cavity or a crack in a piezoelectric material＼[J＼]。 International Journal of Solids and Structures， 1998， 35（17）： 2121-2149。