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

（吉林大學(xué) 機械科學(xué)與工程學(xué)院，吉林 長春 130025）

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

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

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

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

1 機電虛擬裂紋閉合法

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

圖1 機電虛擬裂紋閉合法示意圖Fig．1 Electromechanical virtual crack closure technique

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

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

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

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

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

式中：u1和v1，u2和v2分別為節(jié)點1和節(jié)點2在整體坐標系（X，Y）下位移分量；φ1和φ2分別為節(jié)點1和節(jié)點2在整體坐標系（X，Y）下電勢；Kx和Ky分別為力場下X和Y方向的彈簧剛度；Kφ為電場的彈簧剛度．

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

式中：（x1，y1）和（x5，y5）分別為節(jié)點1和節(jié)點5在整體坐標系（X，Y）下的坐標．如果Δa在每個增量步中都通過位移來更新，裂紋方向隨之更新，在大變形分析問題時很有用．

圖2 啞節(jié)點斷裂壓電單元Fig．2 Fracture of piezoelectric element with dummy nodes

張開位移和電勢差為：

能量釋放率可由式（12）和式（13）代入式（5）－式（7）計算得到．

3 數(shù)值算例

如圖3所示，一含中心裂紋的壓電體，裂紋長度為2a，壓電材料的極化方向為P，邊長2l＝60cm，受均勻拉伸σ∞＝1×105Pa和電位移D∞＝7．4×10－5C／m2的作用，采用PZT－4，P7和PZT－5H3種材料進行數(shù)值模擬，材料屬性見表1，對于該裂紋問題，能量釋放率的理論解為［16］．

其中E∞與D∞關(guān)系為：

圖3 含中心裂紋壓電體模型Fig．3 Piezoelectric model with a centre crack

由于結(jié)構(gòu)和載荷的對稱性，取1／4結(jié)構(gòu)進行求解，頂部施加相應(yīng)的應(yīng)力和電位移，約束左端所有節(jié)點的水平方向位移和底部裂尖以右的節(jié)點電勢（即底部電勢為零）和豎直方向位移，為驗證機電虛擬裂紋閉合法（EMVCCT）的可靠性，當裂紋長度2a＝2 cm 時，將結(jié)構(gòu)離散為I型（單元：15×15），Ⅱ型（單元：30×30）和Ⅲ型（單元：60×60）3種均勻分布網(wǎng)格形式，該結(jié)構(gòu)為I型裂紋，僅需考慮含裂紋壓電體的GI和GD，表2為EMVCCT 計算得到的GI，GD和理論解．由表2可以看出，EMVCCT 在3種網(wǎng)格離散形式、3種材料下均得到了精度較高的GI和GD，與理論解誤差最大不超過3．04%，3種網(wǎng)格所得精度基本一致，可見，該方法不僅具有較高的精度，且對網(wǎng)格的尺寸大小不敏感．

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

表1 材料屬性Tab．1 Material constants

表2 不同網(wǎng)格劃分下應(yīng)變能釋放率Tab．2 Strain energy release rate under different mesh dividing methods

表3 不同裂紋長度下應(yīng)變能釋放率Tab．3 Strain energy release rate under different crack lengths （10－3 N·m－1）

4 結(jié) 論

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

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

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

［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 strain－energy 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 web－stiffener of ship structural details by an objec－ted－oriented 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é)報：自然科學(xué)版，2014，41（8）：13－17．

ZHOU Li－ming，MENG Guang－wei，WANG Hui，etal．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 法的動態(tài)斷裂力學(xué)問題［J］．湖南大學(xué)學(xué)報：自然科學(xué)版，2012，39（11）：41－45．

LONG Shu－yao，ZHANG Guo－h(huán)u．An analysis of the dynamic fracture problem by the meshless local Petrov－Galerkin method［J］．Journal of Hunan University：Natural Sciences，2012，39（11）：41－45．（In Chinese）

［13］平學(xué)成，陳夢成，謝基龍，等．基于新型裂尖雜交元的壓電材料斷裂力學(xué)［J］．力學(xué)學(xué)報，2006，38（3）：407－413．

PING Xue－cheng，CHEN Meng－cheng，XIE Ji－long，etal．Fracture mechanics researches on piezoelectric materials based on a novel crack－tip 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，etal．Cell－based smoothed finite element method－virtual 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 electro－elastic analysis of a cavity or a crack in a piezoelectric material［J］．International Journal of Solids and Structures，1998，35（17）：2121－2149．