基于ABAQUS平臺的機電虛擬裂紋閉合法＊

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

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

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

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

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

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

1 機電虛擬裂紋閉合法

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

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

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

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

2 啞節點斷裂壓電單元

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

裂紋尖端處的節點力和節點電荷為：

式中：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之間的距離：

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

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

張開位移和電勢差為：

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

3 數值算例

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

其中E∞與D∞關系為：

圖3 含中心裂紋壓電體模型Fig．3 Piezoelectric 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 具有表達式簡單，容易編程，只需一步有限元分析，計算精度高的優點．

表1 材料屬性Tab．1 Material constants

表2 不同網格劃分下應變能釋放率Tab．2 Strain energy release rate under different mesh dividing methods

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

4 結 論

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

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 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］周立明，孟廣偉，王暉，等．基于光滑有限元的含裂紋復合材料的虛擬裂紋閉合法［J］．湖南大學學報：自然科學版，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 法的動態斷裂力學問題［J］．湖南大學學報：自然科學版，2012，39（11）：41－45．

LONG Shu－yao，ZHANG Guo－hu．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］平學成，陳夢成，謝基龍，等．基于新型裂尖雜交元的壓電材料斷裂力學［J］．力學學報，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．