基于ABAQUS平臺(tái)的機(jī)電虛擬裂紋閉合法

周立明+孟廣偉+李鵬+李鋒+李宵琳

摘要：為提高求解含裂紋壓電體的能量釋放率的計(jì)算效率和精度以及減少程序編寫與調(diào)試的工作量，針對(duì)四節(jié)點(diǎn)平面壓電單元提出了機(jī)電虛擬裂紋閉合法計(jì)算公式。基于通用商業(yè)有限元軟件ABAQUS，開發(fā)了啞節(jié)點(diǎn)斷裂壓電單元，編寫用戶自定義子程序UEL，該單元可獨(dú)立求解總能量釋放率分量，對(duì)不同材料和裂紋長度的含裂紋壓電體的能量釋放率分量進(jìn)行了求解，與理論解作了對(duì)比，討論了不同形式網(wǎng)格對(duì)求解精度的影響。數(shù)值算例結(jié)果表明，啞節(jié)點(diǎn)斷裂壓電單元具有精度高、簡單方便和對(duì)網(wǎng)格尺寸不敏感等優(yōu)點(diǎn)。

關(guān)鍵詞：裂紋；機(jī)電虛擬裂紋閉合法；能量釋放率；啞節(jié)點(diǎn)斷裂壓電單元； ABAQUS

中圖分類號(hào)：TB115 文獻(xiàn)標(biāo)識(shí)碼：A

越來越多的壓電器件在苛刻環(huán)境條件下服役，裂紋是導(dǎo)致構(gòu)件失效的主要因素，對(duì)構(gòu)件進(jìn)行斷裂分析的第一步便是斷裂參數(shù)的求解。目前計(jì)算斷裂參數(shù)的方法＼[1-2＼]有外推法、J積分、擴(kuò)展有限元法、虛擬裂紋擴(kuò)展法和虛擬裂紋閉合法。外推法要求裂尖處具有特別細(xì)的網(wǎng)格或采用奇異元；J積分表達(dá)式繁瑣，不容易被工程師廣泛采用；擴(kuò)展有限元法裂尖處單元需采用裂尖漸近位移場(chǎng)函數(shù)進(jìn)行加強(qiáng)；虛擬裂紋擴(kuò)展法需要兩次有限元分析；虛擬裂紋閉合法對(duì)網(wǎng)格尺寸不敏感，表達(dá)式簡單，容易編程，只需一步有限元分析，計(jì)算精度高，在解決實(shí)際工程問題中發(fā)揮著重要作用。

Rybicki和Kanninen＼[3＼]于1977年提出虛擬裂紋閉合技術(shù)，將該方法應(yīng)用于求解含裂紋結(jié)構(gòu)的應(yīng)變能釋放率。Raju＼[4＼]和Xie等＼[5＼]對(duì)虛擬裂紋閉合法進(jìn)行了數(shù)學(xué)解釋。Xie等＼[6-8＼]提出了啞節(jié)點(diǎn)斷裂單元，為虛擬裂紋閉合法的發(fā)展作出了貢獻(xiàn)。虛擬裂紋閉合法已應(yīng)用于復(fù)合材料、功能梯度材料、加強(qiáng)結(jié)構(gòu)的斷裂分析與評(píng)估、巖土材料等領(lǐng)域＼[9-10＼]，但機(jī)電耦合場(chǎng)下求解結(jié)構(gòu)斷裂參數(shù)的虛擬裂紋閉合法還未見報(bào)道。

工程結(jié)構(gòu)由于幾何形狀、材料屬性和加載方式的復(fù)雜性，不得不依賴于數(shù)值方法。隨著科學(xué)技術(shù)的發(fā)展，用數(shù)值方法求解斷裂參數(shù)變得切實(shí)可行，很多計(jì)算方法被用來求解斷裂參數(shù)，例如有限差分法＼[11＼]、無網(wǎng)格方法＼[12＼]、雜交元＼[13＼]、邊界元＼[14＼]和光滑有限元＼[15＼]等，但由于缺少商業(yè)軟件的支持，這些數(shù)值方法的工程實(shí)際應(yīng)用相對(duì)缺乏。選用有限元軟件ABAQUS為平臺(tái)，可直接從軟件計(jì)算的結(jié)果中提取相關(guān)信息，通過編寫用戶自定義單元（UEL）子程序，實(shí)現(xiàn)含裂紋壓電材料斷裂參數(shù)的計(jì)算，可極大地減少程序編寫和調(diào)試的工作量，程序一旦得到驗(yàn)證，很方便應(yīng)用到工程實(shí)際的斷裂分析中，提高程序的通用性和計(jì)算效率。

湖南大學(xué)學(xué)報(bào)（自然科學(xué)版）2015年

第10期周立明等：基于ABAQUS平臺(tái)的機(jī)電虛擬裂紋閉合法

在含裂紋壓電材料中，裂紋擴(kuò)展單位長度所需要的能量稱為總能量釋放率，包含機(jī)械能釋放率GM和電能釋放率GD，GM為裂紋尖端前方的應(yīng)力在裂紋改變位移上所做的功，GD為裂紋尖端前方的電位移在裂紋變化后裂紋面電勢(shì)差上所做的功。本文針對(duì)四節(jié)點(diǎn)平面壓電單元提出了機(jī)電虛擬裂紋閉合法計(jì)算公式。以通用商業(yè)有限元軟件ABAQUS為平臺(tái)，開發(fā)了啞節(jié)點(diǎn)斷裂壓電單元，編寫用戶自定義子程序UEL，該單元可獨(dú)立求解總能量釋放率分量，對(duì)不同材料和裂紋長度的含裂紋壓電體的總能量釋放率分量進(jìn)行了求解，并與理論解作了對(duì)比，討論了不同形式網(wǎng)格對(duì)求解精度的影響。

1機(jī)電虛擬裂紋閉合法

針對(duì)求解含裂紋壓電材料的總能量釋放率的需要，提出了針對(duì)四節(jié)點(diǎn)平面壓電單元的機(jī)電虛擬裂紋閉合法計(jì)算公式?？紤]一含裂紋壓電體，如圖1所示，裂紋長度為a，厚度為B，總能釋放率G為產(chǎn)生面積為ΔA的新裂紋面所需要的能量，Δa為裂紋擴(kuò)展量，于是

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分別為原始裂紋尖端處的法向應(yīng)力、切向應(yīng)力和法向電位移；Δv（2），Δu（2）和Δφ（2）分別為裂紋虛擬擴(kuò)展到a+Δa時(shí)裂紋面上的張開位移、相對(duì)滑動(dòng)位移和電勢(shì)差；GI和GII為機(jī)械能釋放率分量；GD為電能釋放率。

圖1機(jī)電虛擬裂紋閉合法示意圖

Fig。1Electromechanical virtual crack closure technique

基于勢(shì)能的改變與將裂紋閉合一個(gè)擴(kuò)展增量所需的功等效，提出機(jī)電虛擬裂紋閉合法。將虛擬裂紋閉合算法在其提出的假設(shè)基礎(chǔ)上進(jìn)行了橫向擴(kuò)展，即在虛擬裂紋擴(kuò)展過程中計(jì)入電勢(shì)以及位移的作用，相對(duì)應(yīng)的壓電單元，將電勢(shì)作為一個(gè)“位移”分量進(jìn)行考慮，式（2）-式（4）可改寫為：

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

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

GD=Δ2BΔa。（7）

式中：x和y，Δ和Δ，，Δ分別為局部坐標(biāo)系（，）下節(jié)點(diǎn)力分量、張開位移分量、節(jié)點(diǎn)電荷和電勢(shì)差。

由式（5）-式（7）可知：1）機(jī)電虛擬裂紋閉合法可分別計(jì)算GI，GII，GD；2）能量釋放率的計(jì)算僅僅包含節(jié)點(diǎn)力與節(jié)點(diǎn)位移、節(jié)點(diǎn)電荷與節(jié)點(diǎn)電勢(shì)差，這些變量可從有限元軟件中輸出；3）避免了應(yīng)力和電位移的積分，容易與有限元分析相結(jié)合。

2 啞節(jié)點(diǎn)斷裂壓電單元

基于ABAQUS平臺(tái)，利用啞節(jié)點(diǎn)斷裂壓電單元來實(shí)現(xiàn)二維線狀含裂紋壓電斷裂力學(xué)問題，通過編寫用戶自定義單元子程序UEL來實(shí)現(xiàn)。啞節(jié)點(diǎn)斷裂壓電單元的定義及其節(jié)點(diǎn)編號(hào)如圖2所示，含5個(gè)節(jié)點(diǎn)。節(jié)點(diǎn)1和節(jié)點(diǎn)2對(duì)應(yīng)于裂紋尖端，在其節(jié)點(diǎn)間放置有特殊剛度的彈簧，節(jié)點(diǎn)3和節(jié)點(diǎn)4位于裂紋尖端的后面，節(jié)點(diǎn)5在裂紋尖端的前面。在ABAQUS中，單元所具有的完整的節(jié)點(diǎn)矢量為：

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

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

裂紋尖端處的節(jié)點(diǎn)力和節(jié)點(diǎn)電荷為：

Fx=Kx（u1

Symbolm@@ u2）， Fy=Ky（v1

Symbolm@@ v2），

Q=K（1

Symbolm@@ 2）。（9）

式中：u1和v1，u2和v2分別為節(jié)點(diǎn)1和節(jié)點(diǎn)2在整體坐標(biāo)系（X，Y）下位移分量；1和2分別為節(jié)點(diǎn)1和節(jié)點(diǎn)2在整體坐標(biāo)系（X，Y）下電勢(shì)；Kx和Ky分別為力場(chǎng)下X和Y方向的彈簧剛度；K為電場(chǎng)的彈簧剛度。

節(jié)點(diǎn)3、節(jié)點(diǎn)4和節(jié)點(diǎn)5被用來從ABAQUS結(jié)果中提取相關(guān)信息，對(duì)單元的剛度矩陣并沒有實(shí)際貢獻(xiàn)，稱為“啞節(jié)點(diǎn)”，該單元為啞節(jié)點(diǎn)斷裂壓電單元。虛擬裂紋擴(kuò)展量是主節(jié)點(diǎn)1和啞節(jié)點(diǎn)5之間的距離：

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

式中：（x1，y1）和（x5，y5）分別為節(jié)點(diǎn)1和節(jié)點(diǎn)5在整體坐標(biāo)系（X，Y）下的坐標(biāo)。如果

SymbolDA@ a在每個(gè)增量步中都通過位移來更新，裂紋方向隨之更新，在大變形分析問題時(shí)很有用。

圖2啞節(jié)點(diǎn)斷裂壓電單元

Fig。2Fracture of piezoelectric element

with dummy nodes

應(yīng)變能釋放率必須在裂紋尖端處的局部坐標(biāo)系（，）下計(jì)算，X軸和軸之間的夾角為：

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

通過簡單的矢量投影關(guān)系，在局部坐標(biāo)系（，）下節(jié)點(diǎn)力和節(jié)點(diǎn)電荷為：

x=Fxcosθ+Fysinθ；

y=-Fxsinθ+Fycosθ；

=Q。（12）

張開位移和電勢(shì)差為：

Δ=Δucosθ+Δvsinθ；

Δ=-Δusinθ+Δvcosθ；

=Δφ。（13）

能量釋放率可由式（12）和式（13）代入式（5）-式（7）計(jì)算得到。

3數(shù)值算例

如圖3所示，一含中心裂紋的壓電體，裂紋長度為2a，壓電材料的極化方向?yàn)镻，邊長2l=60 cm，受均勻拉伸

SymbolsA@

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

SymboleB@ =7。4×10-5C/m2的作用，采用PZT4，P7和PZT5H 3種材料進(jìn)行數(shù)值模擬，材料屬性見表1，對(duì)于該裂紋問題，能量釋放率的理論解為＼[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@ 關(guān)系為：

D

SymboleB@ =c11c33-c13e31c11c33-c213σ

SymboleB@ +

d33+c33e231-2c13e31e33+c11e233c11c33-c213E

SymboleB@ 。（17）

圖3含中心裂紋壓電體模型

Fig。3Piezoelectric model with a centre crack

由于結(jié)構(gòu)和載荷的對(duì)稱性，取1/4結(jié)構(gòu)進(jìn)行求

解，頂部施加相應(yīng)的應(yīng)力和電位移，約束左端所有節(jié)點(diǎn)的水平方向位移和底部裂尖以右的節(jié)點(diǎn)電勢(shì)（即底部電勢(shì)為零）和豎直方向位移，為驗(yàn)證機(jī)電虛擬裂紋閉合法（EMVCCT）的可靠性，當(dāng)裂紋長度2a=2 cm時(shí)，將結(jié)構(gòu)離散為I型（單元：15×15），Ⅱ型（單元：30×30）和Ⅲ型（單元：60×60） 3種均勻分布網(wǎng)格形式，該結(jié)構(gòu)為I型裂紋，僅需考慮含裂紋壓電體的GI和GD，表2為EMVCCT計(jì)算得到的GI，GD和理論解。

由表2可以看出，EMVCCT在3種網(wǎng)格離散形式、3種材料下均得到了精度較高的GI和GD，與理論解誤差最大不超過3。04%，3種網(wǎng)格所得精度基本一致，可見，該方法不僅具有較高的精度，且對(duì)網(wǎng)格的尺寸大小不敏感。

為方便計(jì)算不同裂紋長度下結(jié)構(gòu)的GI和GD，采用Ⅲ型單元離散形式，表3給出了3種材料在不同裂紋長度下得到的GI和GD，并與理論解做了比較，從結(jié)果可以看出EMVCCT的計(jì)算結(jié)果比解析解得到的結(jié)果要小，是由于插值函數(shù)使用了“協(xié)調(diào)和完整的位移函數(shù)”，連續(xù)體離散后剛度會(huì)有所增加，求解值相對(duì)實(shí)際值要小，因此，EMVCCT的計(jì)算結(jié)果比解析解得到的結(jié)果要??；從結(jié)果還可以看出EMVCCT具有較高的精度，進(jìn)一步驗(yàn)證了EMVCCT具有表達(dá)式簡單，容易編程，只需一步有限元分析，計(jì)算精度高的優(yōu)點(diǎn)。

4結(jié)論

本文針對(duì)四節(jié)點(diǎn)平面壓電單元提出了機(jī)電虛擬裂紋閉合法，基于通用商業(yè)有限元軟件ABAQUS平臺(tái)，開發(fā)了啞節(jié)點(diǎn)斷裂壓電單元，編寫用戶自定義子程序UEL，以PZT4和P7、PZT5H的壓電平板的中心裂紋問題為例，求解了不同網(wǎng)格離散形式和裂紋長度下結(jié)構(gòu)的總能量釋放率，并與理論解做了對(duì)比，結(jié)論如下：

1）該方法表達(dá)式簡單，容易編程，只需一步有限元分析，對(duì)網(wǎng)格的尺寸大小不敏感，具有較高的計(jì)算精度。

2）該方法基于有限元軟件ABAQUS可直接從軟件計(jì)算的結(jié)果中提取相關(guān)信息，極大地減少了程序編寫和調(diào)試的工作量。

參考文獻(xiàn)

[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]周立明，孟廣偉，王暉，等。 基于光滑有限元的含裂紋復(fù)合材料的虛擬裂紋閉合法＼[J＼]。 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版，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法的動(dòng)態(tài)斷裂力學(xué)問題＼[J＼]。 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版，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]平學(xué)成， 陳夢(mèng)成， 謝基龍，等。 基于新型裂尖雜交元的壓電材料斷裂力學(xué)＼[J＼]。 力學(xué)學(xué)報(bào)， 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。