雷 鈞，劉和國，楊慶生

（北京工業(yè)大學(xué)機(jī)械工程與應(yīng)用電子技術(shù)學(xué)院，北京 100124）

基于相互作用積分法的壓電裂紋邊界元分析

雷 鈞，劉和國，楊慶生

（北京工業(yè)大學(xué)機(jī)械工程與應(yīng)用電子技術(shù)學(xué)院，北京 100124）

為了獲得更為精確高效的壓電裂紋分析方法，推導(dǎo)了壓電材料中基于相互作用積分計(jì)算裂尖場(chǎng)集中系數(shù)的計(jì)算公式，并成功應(yīng)用于相應(yīng)的邊界元分析軟件.應(yīng)用此軟件，對(duì)PZT-4壓電材料在不同加載條件、不同的極化方向等情況下壓電矩形板中的中心直裂紋以及邊裂紋問題進(jìn)行了分析，驗(yàn)證了相互作用積分法的路徑無關(guān)性，并且計(jì)算出各型裂尖場(chǎng)集中系數(shù)，與相應(yīng)的基于J積分法以及位移插值法獲得的結(jié)果對(duì)比良好，說明該方法是一種高效準(zhǔn)確的壓電裂紋分析方法.

相互作用積分；邊界元法；壓電材料；強(qiáng)度因子；極化方向

壓電材料因其獨(dú)特的壓電效應(yīng)在智能材料與結(jié)構(gòu)中得到了廣泛的應(yīng)用.但常見的壓電陶瓷為脆性材料，斷裂韌性低、缺陷敏感性高，在工作載荷下易發(fā)生斷裂失效.因而對(duì)其進(jìn)行斷裂研究十分重要.

目前針對(duì)壓電材料的斷裂分析包括實(shí)驗(yàn)研究、理論分析以及數(shù)值模擬等.理論研究主要針對(duì)無限大壓電體等特殊情況，如Suo等［1］和Pak［2］對(duì)無限大壓電體絕緣裂紋平面問題給出了解析解.實(shí)驗(yàn)研究則主要針對(duì)靜態(tài)斷裂問題，如Park等［3-4］在1995年先后對(duì)PZT-4緊湊拉伸試件進(jìn)行了斷裂實(shí)驗(yàn)及三點(diǎn)彎實(shí)驗(yàn).數(shù)值模擬方便快捷，具有重要的現(xiàn)實(shí)意義.Enderlin等［5］實(shí)現(xiàn)了壓電J積分在有限元法中的運(yùn)用.Rao等［6］和Yu等［7］實(shí)現(xiàn)了相互作用積分法在有限元法對(duì)壓電材料斷裂研究中的應(yīng)用.Lei等［8-11］實(shí)現(xiàn)了 J積分法［12］、位移法、公式應(yīng)力法在邊界元法對(duì)壓電材料斷裂參數(shù)計(jì)算中的應(yīng)用，數(shù)值驗(yàn)證了J積分法計(jì)算壓電斷裂參數(shù)的高精度.如Lei等［8］文中所述，由J積分計(jì)算各型場(chǎng)集中系數(shù)，有2種基本方法:一種是通過基于裂尖張開位移比率的插值公式獲得；另一種是基于輔助場(chǎng)的解析公式提出的所謂相互作用積分法［13-14］獲得.前者由于需要計(jì)算裂尖張開位移，這顯然依賴于數(shù)值計(jì)算高應(yīng)變梯度裂尖區(qū)域的應(yīng)變分布精度，對(duì)于有限元等區(qū)域型數(shù)值方法而言，會(huì)嚴(yán)重依賴于局部網(wǎng)格劃分.而后者相互作用積分法無需計(jì)算裂紋面參數(shù)，裂尖單元的選取對(duì)計(jì)算結(jié)果影響較小，在計(jì)算精度上更具優(yōu)勢(shì).而與傳統(tǒng)J積分法相比，相互作用積分法保留了J積分路徑無關(guān)性的特點(diǎn)，同時(shí)與J積分法不同的是相互作用積分法可以直接得到各型應(yīng)力強(qiáng)度因子，無需進(jìn)一步分離，計(jì)算更為簡(jiǎn)單快捷.與傳統(tǒng)位移法、公式應(yīng)力法相比，相互作用積分法無需在裂紋面進(jìn)行插值計(jì)算，計(jì)算結(jié)果具有更高的精度.J積分公式在有限元實(shí)現(xiàn)時(shí)，需引入一個(gè)光滑的權(quán)函數(shù)將圍道積分轉(zhuǎn)化為等效的區(qū)域積分，其迭代插值得到的計(jì)算結(jié)果對(duì)單元?jiǎng)澐志哂懈呙舾卸?而邊界元法可直接進(jìn)行J積分中的圍道積分，無需插值迭代，相對(duì)簡(jiǎn)單高效.而據(jù)作者所知，目前，還未見在邊界元法中采用相互作用積分法計(jì)算壓電斷裂參數(shù)的相關(guān)文獻(xiàn).因此，本文將壓電材料中基于相互作用積分公式，應(yīng)用于相應(yīng)的邊界元分析軟件，實(shí)現(xiàn)對(duì)壓電材料中的裂尖斷裂參數(shù)的計(jì)算，并與J積分、位移法等計(jì)算結(jié)果進(jìn)行對(duì)比.

1　壓電方程

令σ、u、D分別表示應(yīng)力張量、彈性張量和電位移張量.在準(zhǔn)靜電場(chǎng)假設(shè)以及無體積力和體電荷情況下，二維線性均勻、各向異性壓電體的平衡方程如下:

式中Δ表示哈密頓算子.

令ε、E分別表示應(yīng)變張量和電場(chǎng)張量，它們分別為位移u和電勢(shì)φ的梯度

壓電體的本構(gòu)方程可表示為

式中:C、e、κ分別為4階彈性張量、3階壓電張量和2階介電張量，三者滿足對(duì)稱性

2　邊界積分方程

對(duì)于二維壓電裂紋體，其外部邊界記為?Ω，裂紋面ΓC由上裂紋面Γ+C與下裂紋面Γ-C組成，即ΓC=Γ-C+Γ+C.引入廣義位移張量u=［ui，φ］T，應(yīng)力σ=［σij，Di］T，應(yīng)變 ε=［εij，-Ei］T，面力 p= σiJn，并定義廣義張開位移為ΔuJ=uJ（x∈Γ+C）-uJ（x∈Γ-C）.

假設(shè)壓電裂紋滿足Deeg［15］提出的絕緣型裂紋面電邊界條件，借助Lei等［8］文中的對(duì)偶邊界積分方程法，分別對(duì)外邊界?Ω與裂紋面Γ+C應(yīng)用位移及面力邊界積分方程式中:x和y分別表示源點(diǎn)和觀察點(diǎn)；cIJ=δIJ/2是邊界平滑參數(shù)；uGIJ和pGIJ分別為二維壓電體靜態(tài)位移與面力基本解；dGIJ和 sGIJ表示其高階基本解，具體表達(dá)式可參考文獻(xiàn)［16］.

3　相互作用積分

裂紋尖端的局部應(yīng)變能釋放率可由與積分路徑無關(guān)的J積分公式計(jì)算獲得，二維均勻壓電介質(zhì)的J積分表達(dá)式由Pak［12］給出:

式中:Γ為圍繞裂尖的任意閉合積分路徑；nj表示閉合積分路徑的單位外法向量分量，W=σiJεiJ/2為線性壓電的電焓密度.

考慮壓電斷裂模型2個(gè)獨(dú)立的平衡狀態(tài):一個(gè)是基于實(shí)際加載條件下的真實(shí)場(chǎng)狀態(tài)；另外一個(gè)是適用于任一斷裂模型裂尖場(chǎng)的輔助場(chǎng)狀態(tài).定義在真實(shí)載荷與輔助場(chǎng)聯(lián)合作用下的J積分記為JTol，表達(dá)式為

其中上標(biāo)“aux”表示輔助場(chǎng)的力學(xué)與電學(xué)變量，

I表示相互作用積分，應(yīng)用等式

相互作用積分I可表示為

Rao和Kuan［6］定義二維Griffith裂紋尖端場(chǎng)作為輔助場(chǎng)變量，具體表達(dá)式為

式中:裂尖場(chǎng)的強(qiáng)度因子K=［K?，KⅠ，KD］T；fNiJ（θ）和dNI（θ）表示均勻壓電體的標(biāo)準(zhǔn)角函數(shù)，具體表達(dá)式見文獻(xiàn)［6］.

在數(shù)值計(jì)算時(shí)，可取圍道積分路徑，如圖1所示.

進(jìn)一步將圍道離散成N個(gè)二次連續(xù)單元，每個(gè)單元取的高斯點(diǎn)的個(gè)數(shù)為k，則式（11）離散為

式中:L為雅克比矩陣；wm表示高斯點(diǎn)的權(quán)重.

與J積分類似，壓電材料的相互作用積分可表示為

式中:Y為與材料參數(shù)的有關(guān)的Irwin矩陣，特殊地，如果取?型裂紋的輔助場(chǎng)，即Kaux

?=1，Kaux

Ⅰ=KauxD= 0，式（15）滿足

同理，對(duì)于其他類型裂紋，可得

綜上所述，廣義應(yīng)力強(qiáng)度因子可由求解下列方程得到:

4　數(shù)值算例

以下的算例中以常見的PZT-4壓電材料作為研究對(duì)象，其材料參數(shù)為

裂尖場(chǎng)強(qiáng)度因子的量綱一計(jì)算公式為

4.1矩形壓電板中心直裂紋強(qiáng)度因子計(jì)算

第1個(gè)算例取如圖2所示含中心直裂紋的PZT-4矩形壓電板，具體尺寸為:2W=1.0 m，2L= 2.0 m，2a=0.2 m，在x2方向施加載荷（σ022，D02），令加載參數(shù)為λ=c22D02/（e22σ022），取圖3所示3個(gè)圍道積分路徑（以右裂尖為中心，邊長(zhǎng)分別為0.18、0.2、0.22 m的正方形圍道路徑1、路徑2、路徑3，在表1中記為P1、P2、P3），分別通過相互作用積分計(jì)算不同加載條件下的強(qiáng)度因子，并將結(jié)果與J積分及位移法結(jié)果進(jìn)行對(duì)比，如表1所示.相互作用積分結(jié)果與J積分結(jié)果最大誤差小于0.25%，充分證明了相互作用積分的高效準(zhǔn)確性與路徑無關(guān)性.

對(duì)上述含中心裂紋的PZT-4壓電板，取路徑2為計(jì)算路徑，分別對(duì)每個(gè)二次單元取不同的高斯點(diǎn)數(shù)，力電強(qiáng)度因子變化情況如圖4所示.由圖4可以看出，在不同加載條件下，高斯點(diǎn)的選取對(duì)壓電板不同加載條件下力電斷裂參數(shù)的計(jì)算結(jié)果基本沒有影響.

4.2矩形壓電板邊裂紋強(qiáng)度因子計(jì)算

圖5所示含邊裂紋的矩形壓電板，極化方向?yàn)檠刂鴛2軸方向，具體尺寸為:2W=5 cm，2L=4 cm，2a=1.2 cm，由于算例4.2中已經(jīng)論證了相互作用積分的路徑無關(guān)性，為減少計(jì)算量，取圖5所示2條圍道作為積分路徑（以右裂尖為中心，邊長(zhǎng)分別為1.2 cm和1.6 cm的正方形圍道P1、P2），每個(gè)單元取20個(gè)高斯點(diǎn)，分別計(jì)算在不同加載條件下含邊裂紋PZT-4壓電板的斷裂參數(shù)，并將相關(guān)結(jié)果與位移法、J積分法結(jié)果進(jìn)行對(duì)比，相互作用積分結(jié)果與J積分結(jié)果最大誤差小于0.6%，充分證明了相互作用積分的高效準(zhǔn)確性與路徑無關(guān)性，如表2所示.

表1　中心直裂紋λ=-7.6,0,7.6加載斷裂參數(shù)Table 1 Fracture parameter for the center crack under λ=-7.6,0,7.6 loading

表2　邊裂紋λ=-7.6、0、7.6加載斷裂參數(shù)Table 2 Fracture parameter for the edge crack under λ=-7.6,0,7.6 loading

4.3極化角度對(duì)斷裂韌性的影響

對(duì)上述含中心裂紋的PZT-4壓電板，取極化方向垂直于裂紋面方向?yàn)?°，于極化方向?yàn)?0°，分別計(jì)算極化角度為0°～90°PZT-4壓電板在力電耦合（σ022=1 GPa，D02=1 C/m2）作用下的斷裂參數(shù)力電強(qiáng)度因子變化情況如表3所示.為減少計(jì)算量，相互作用積分結(jié)果只與位移法結(jié)果進(jìn)行了對(duì)比.

由表3中數(shù)據(jù)可知，2種方法對(duì)比良好，最大誤差不超過0.2%.2種計(jì)算方法均表明在相同加載條件下:PZT-4壓電板在極化0°時(shí)電位移強(qiáng)度因子最大，而應(yīng)力位移強(qiáng)度因子隨著極化角度變化基本保持不變.

表3　極化角度對(duì)PZT-4壓電板斷裂參數(shù)的影響Table 3 Influence of polarization angle on the fracture parameter

5　結(jié)論

1）在計(jì)算壓電材料斷裂參數(shù)時(shí)，本文提出的基于邊界元法的相互作用積分是精確高效的.

2）采用相互作用積分計(jì)算含中心裂紋及邊裂紋PZT-4壓電板在不同加載情況下的斷裂參數(shù)時(shí)與J積分、位移法對(duì)比良好，且計(jì)算結(jié)果與圍道積分路徑無關(guān).

3）相互作用積分采取了圍道積分進(jìn)行計(jì)算，圍道積分上每個(gè)單元選取的高斯點(diǎn)的個(gè)數(shù)對(duì)不同加載情況下壓電板斷裂參數(shù)計(jì)算結(jié)果影響不大.

4）相互作用積分法與位移法均表明PZT-4壓電板在極化0°時(shí)電位移強(qiáng)度因子最大，而應(yīng)力位移強(qiáng)度因子隨著極化角度變化基本保持不變.

［1］SUO Z G，KUO C M，BAMETT D M，et al.Fracture mechanics for piezoelectric ceramics［J］.Journal of the Mechanics and Physics of Solids，1992，40（4）:739-765.

［2］PAK Y E.Linear electroelatic fracture-mechanics of piezoelectricmaterials［J］.InternationalJournalofFracture，1992，54（1）:79-100.

［3］PARK S B，SUN C T.Fracture criteria for piezoelectric ceramics［J］.Journal of the American Ceramic Society，1995，78（6）:1475-1480.

［4］PARK S B，SUN C T.Effect of electric-field on fracture of piezoelectricceramics［J］.InternationalJournalof Fracture，1995，70（3）:203-216.

［5］ENDERLIN M，RICOEUR A，KUNA M.Finite element techniques for dynamic crack analysis in piezoelectrics ［J］.International Journal of Fracture，2005，134（3）: 191-208.

［6］RAO B N，KUNA M.Interaction integrals for fracture analysis of functionally graded piezoelectric materials［J］. International Journal of Solids and Structures，2008，45 （20）:5237-5257.

［7］YU H J，WU L Z，GUO L C，et al.A domainindependent interaction integral for fracture analysis of nonhomogeneous piezoelectric materials［J］.International Journal of Solids and Structures，2012，49（23）:3301-3315.

［8］LEI J，WANG H Y，ZHANG C Z，et al.Comparison of several BEM-based approaches in evaluating crack-tip field intensity factors in piezoelectric material［J］.International Journal of Fracture，2014，189（1）:111-120.

［9］LEI J，ZHANG C Z，F(xiàn)ELIPE G S.BEM analysis of electrically limited permeable cracks considering Coulomb tractions inpiezoelectricmaterials［J］.Engineering Analysis with Boundary Elements，2015，54:28-38.`

［10］LEI J，WANG H Y，ZHANG C Z，et al.Crack-tip amplificationandshieldingbymicro-cracksin piezoelectric solids-Part I:Static case［J］.Engineering Analysis with Boundary Elements，2013，37（12）:1585-1593.

［11］LEI J，YANG Q S，ZHANG C Z.Transient response of a semi-permeablecrackbetweendissimilaranisotropic piezoelectriclayersbytime-domainBEM［J］. Engineering Analysis with Boundary Elements，2013，37 （10）:1205-1211.

［12］PAK Y E.Crack extension force in a piezoelectric material［J］.Journal of Applied Mechanics，1990，57 （3）:647-653.

［13］RAO B N，RAHMAN S.Mesh-free analysis of cracks in isotropic functionally graded materials［J］.Engineering Fracture Mechanics，2003，70（1）:1-27.

［14］RAO B N，RAHMAN S.An interaction integral method for analysis of cracks in orthotropic functionally graded materials［J］.Computational Mechanics 2003，32（1）: 40-51.

［15］DEEG W F J.The analysis of dislocation，crack，and inclusion problem in piezoelectric solids［D］.Stanford: Stanford University，1980.

［16］FELIPE G S，ZHANG C Z，ANDRES S.2-D transient dynamic analysis of cracked piezoelectric solids by a time-domain BEM［J］.Computer Methods in Applied Mechanics and Engineering，2008，197（33）:3108-3121.

（責(zé)任編輯 楊開英）

Analysis of Boundary Element for Piezoelectric Cracks Using an Interaction Integral Method

LEI Jun，LIU Heguo，YANG Qingsheng
（College of Mechanical Engineering and Applied Electronics Technology，Beijing University of Technology，Beijing 100124，China）

To get a more effective and accurate method to study the fracture behavior of the piezoelectric materials，a computational formula based on the interaction integral method to calculate the crack-tip field intensity factors was obtained in this paper.This method was successfully applied to the boundary element analysis software.With this software，the fracture behavior of PZT-4 made piezoelectric board's central crack and edge crack under different loading conditions was analyzed in different polarization directions. The path-independent characteristics of the interaction integral method was verified.The crack-tip field intensity factors obtained were compared well with the results of both the J integral method and the displacement interpolating method.This proves that this method is effective and accurate.

interaction integral；boundary element method；piezoelectric material；intensify factors；polarization direction

U 461；TP308

A

0254-0037（2016）02-0161-06

10.11936/bjutxb2015030009

2015-03-04

國家自然科學(xué)基金資助項(xiàng)目（11472021）

雷 鈞（1976—），男，副教授，主要從事新型智能材料力學(xué)性能評(píng)價(jià)、數(shù)值計(jì)算技術(shù)方面的研究.E-mail:leijun@ bjut.edu.cn