雙材料界面裂紋復(fù)應(yīng)力強(qiáng)度因子的正則化邊界元法1)

谷 巖 張耀明

?(青島大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,山東青島 266071)

?(山東理工大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,山東淄博 255049)

引言

雙材料比較容易在結(jié)合面附近萌生裂縫,并沿著界面發(fā)生擴(kuò)展直至結(jié)構(gòu)破壞.因此,雙材料的破壞不僅表現(xiàn)為單一材料內(nèi)部的裂紋問題,往往更多的是界面的破壞、開裂或脫黏[1].隨著計(jì)算力學(xué)學(xué)科的不斷發(fā)展,用于解決斷裂力學(xué)問題的數(shù)值計(jì)算方法不斷涌現(xiàn),從早期的有限差分法[2]、有限元法[3]、邊界元法[4-9]到現(xiàn)在的無網(wǎng)格法[10-11]、數(shù)值流形法[12]、非連續(xù)變形分析法[13]、分子動力學(xué)法[14]等,它們成為推動斷裂力學(xué)研究不斷發(fā)展的重要工具.

采用數(shù)值方法進(jìn)行斷裂力學(xué)分析時,裂紋尖端奇異區(qū)域處理的好壞直接關(guān)系到最終斷裂力學(xué)參數(shù)(如應(yīng)力強(qiáng)度因子、能量釋放率等)的求解精度.對于傳統(tǒng)單一材料,材料內(nèi)部裂紋尖端位移和應(yīng)力場具有經(jīng)典的平方根(r1/2)和負(fù)平方根(r?1/2)漸近性,相應(yīng)的理論研究和數(shù)值求解技術(shù)已發(fā)展的非常成熟.比如對裂尖奇異應(yīng)力場的求解,已發(fā)展了包括“1/4 單元法”、“數(shù)值外插法”、“虛擬裂紋閉合法”和“J 積分/相互作用積分”等多種有效的求解方法.然而,與傳統(tǒng)單一材料不同,雙材料界面裂紋漸近位移和應(yīng)力場表現(xiàn)出劇烈的振蕩特性,即[15-16]

其中i 為復(fù)數(shù)單位,ε 為雙材料參數(shù)(或振蕩因子).導(dǎo)致許多用于表征經(jīng)典的平方根和負(fù)平方根物理場漸近性的傳統(tǒng)方法失效[17-18].雖然界面裂紋復(fù)應(yīng)力強(qiáng)度因子(K1+iK2)可采用各種守恒積分(如J 積分)或位移外插法等進(jìn)行間接求解,但計(jì)算精度較低,也難以直接獲得裂尖近場區(qū)域物理量的精確分布[19].因此,雖然目前計(jì)算斷裂力學(xué)發(fā)展的比較成熟,但其在雙材料界面裂紋分析方面的研究依然比較欠缺,從基礎(chǔ)研究的角度來考慮還需進(jìn)行大量的研究.

此外,隨著實(shí)際工程需求和表面涂層工藝的逐漸成熟,大量雙材料的相對厚度(表面特征長度和厚度的尺寸比)越來越小,特別是各種微納米復(fù)合材料的相繼出現(xiàn),給精確有效的數(shù)值模擬帶來了巨大的挑戰(zhàn).比如,為了避免出現(xiàn)畸形網(wǎng)格,有限元法需按照涂層的厚度劃分單元.可這樣做必然會導(dǎo)致百萬甚至幾百萬個子單元,計(jì)算工作量劇增,同時解的收斂性也難以保證[20].相對于有限元法,邊界元僅需在邊界離散,可以從容處理邊界單元的疏密過渡,很大程度上避免了有限元法中網(wǎng)格畸變帶來的困難,這一點(diǎn)對薄體與涂層結(jié)構(gòu)問題的數(shù)值模擬尤為重要.近年來,邊界元法在薄體結(jié)構(gòu)問題中的應(yīng)用正在蓬勃發(fā)展,被認(rèn)為是其比有限元法更具優(yōu)勢的一個新領(lǐng)域,而且在表面涂層、壓電薄膜等工程領(lǐng)域展現(xiàn)出了廣闊的應(yīng)用前景[21-27].

論文以邊界元法為基本工具,通過引入一種能精確表征復(fù)合材料界面裂紋振蕩特性的“特殊裂尖單元”[28],實(shí)現(xiàn)了雙材料界面裂紋復(fù)應(yīng)力強(qiáng)度因子的精確求解.該方法無需裂尖區(qū)域高密度的網(wǎng)格剖分,可顯著提高裂尖近場力學(xué)參量和復(fù)應(yīng)力強(qiáng)度因子的求解精度和數(shù)值穩(wěn)定性.在此基礎(chǔ)上,結(jié)合邊界元法中計(jì)算近奇異積分的一種通用正則化算法[29-30],成功求解了超薄雙材料界面裂紋問題.數(shù)值算例表明,所提算法穩(wěn)定,效率高,在不增加計(jì)算量的前提下,顯著提高了裂尖近場力學(xué)參量和斷裂力學(xué)參數(shù)的求解精度和數(shù)值穩(wěn)定性.

1 界面裂紋漸進(jìn)位移和應(yīng)力場的振蕩特性

圖1 中,考慮由兩種不同材料組成的雙材料板,在板的結(jié)合面存在一條界面裂紋.記上層和下層材料分別為區(qū)域?1和?2,相應(yīng)的彈性模量和泊松比分別為μ1,ν1和μ2,ν2.

圖1 雙材料界面裂紋問題Fig.1 An interface crack in a dissimilar material

早在20 世紀(jì)五六十年代,Williams[16],Erdogan[15]和England[17]等就發(fā)現(xiàn)雙材料界面裂紋漸近位移和應(yīng)力場具有形如式(1)的振蕩特性.其中r為計(jì)算點(diǎn)到裂尖的距離,雙材料參數(shù)ε 定義如下

對于平面應(yīng)變問題,有κi=3 ?4νi.對于平面應(yīng)力問題,有κi=(3 ?νi)/(1+νi).裂紋尖端附近的位移和應(yīng)力場可表示為[28,31]

K為復(fù)應(yīng)力強(qiáng)度因子,ci=(1+κi)/μi為材料常數(shù).由式(3)或式(4)可知,復(fù)應(yīng)力強(qiáng)度因子K1+iK2的值可通過裂尖張開位移或應(yīng)力進(jìn)行數(shù)值求解.

傳統(tǒng)有限元方法在處理裂紋問題時,需嚴(yán)格按照材料界面剖分網(wǎng)格,且裂尖區(qū)域需要采用高密度網(wǎng)格以捕捉應(yīng)力場的奇異性,前處理工作比較繁瑣.此外,需要指出,界面裂紋漸近位移和應(yīng)力場具有r1/2+iε和r?1/2+iε的振蕩特性,許多用于表征經(jīng)典的平方根和負(fù)平方根物理場漸近性的傳統(tǒng)方法失效,這也給雙材料界面裂紋精確有效的數(shù)值模擬帶來了不小的挑戰(zhàn).

2 求解界面裂紋問題的邊界元方法

采用多域邊界元法對圖1 所示的雙材料界面裂紋問題進(jìn)行數(shù)值求解.對區(qū)域?1和?2,可分別建立如下邊界積分方程[28]

其中P和Q代表計(jì)算點(diǎn)和邊界積分點(diǎn),uj和tj(j=1,2)代表邊界位移和面力,Uij(P,Q)和Tij(P,Q)為位移和面力基本解.對于光滑邊界,Ci j(P)=δi j/2.為了保證足夠的計(jì)算精度,采用三點(diǎn)二次非連續(xù)單元近似(圖2).邊界元法中的單元近似分為對邊界形狀的“幾何近似”和對邊界物理參量的“物理近似”.對幾何邊界的近似可表示為

圖2 二次非連續(xù)單元Fig.2 Discontinuous quadratic element

將邊界積分方程(6)分別應(yīng)用于區(qū)域?1和?2,可構(gòu)造如下線性代數(shù)方程組

其中H和G代表系數(shù)矩陣,上標(biāo)“1”和“2”分別表示區(qū)域?1和?2,下標(biāo)“I”代表雙材料結(jié)合面(界面)處對應(yīng)的系數(shù)矩陣或物理量.根據(jù)界面位移和面力的連續(xù)條件

方程組(12)和(13)可以耦合為

通過對方程(15) 的求解,可得到邊界點(diǎn)處位移和面力的數(shù)值解.將裂紋尖端的張開位移或面力代入式(3) 或式(4),便可進(jìn)一步求得復(fù)應(yīng)力強(qiáng)度因子的值.

需要指出,界面裂紋漸近位移和應(yīng)力場具有r1/2+iε和r?1/2+iε的振蕩特性.因此,在采用以上分域法的同時,需要對裂尖單元(包含裂紋尖端的左右兩個單元,見圖3)進(jìn)行特殊處理,這就是接下來要引入的“特殊裂尖單元”法.

圖3 界面裂紋裂尖單元Fig.3 Elements near the crack-tip

3 雙材料界面裂紋“特殊裂尖單元”

引入一種含有復(fù)振蕩因子的新型“特殊裂尖單元”[28],可精確表征裂紋尖端漸近位移和應(yīng)力場的振蕩特性.提高界面裂紋復(fù)應(yīng)力強(qiáng)度因子的數(shù)值求解精度.下文簡要介紹方法的數(shù)學(xué)原理.

對于幾何邊界的近似,裂尖單元采用如式(7)和式(8) 相同的二次單元.但對于位移和面力的近似(物理量的近似),裂尖單元需要做如下特殊處理.如式(3) 所示,為了準(zhǔn)確表征裂尖位移場的振蕩特性,裂尖位移場在數(shù)學(xué)上應(yīng)該具有以下形式

其中d1代表剛體位移,d2r1/2+iε為主要項(xiàng),d3r3/2+iε為距離函數(shù)r的高階無窮小量(對計(jì)算精度的影響可以忽略).

數(shù)值上,通過式(9) 插值得到的裂尖位移場應(yīng)具有如式(16)的漸進(jìn)性.當(dāng)裂紋尖端在單元中的無因次坐標(biāo)為ξ=?1 時,即裂尖位于單元左端點(diǎn)時(圖3),式(9) 中的位移插值形函數(shù)應(yīng)表示為以下形式

4 薄體結(jié)構(gòu)邊界元法中近奇異積分的計(jì)算

薄體結(jié)構(gòu)問題的數(shù)值求解是邊界元法的難點(diǎn)之一,其實(shí)質(zhì)是近奇異積分的精確計(jì)算.原因是,受薄體結(jié)構(gòu)特殊幾何構(gòu)造的影響,邊界節(jié)點(diǎn)通常和某些積分單元十分地接近,導(dǎo)致邊界積分方程(6)的積分表現(xiàn)出不同程度的近奇異性,常規(guī)高斯積分結(jié)果失效.本文采用一種非線性變量替換法對近奇異積分進(jìn)行精確求解.該方法基于盡可能拉近運(yùn)算因子間數(shù)量級或變化尺度的思想,可有效改善近奇異核函數(shù)的振蕩特性,使近奇異積分的計(jì)算結(jié)果得到大幅提高.

二維彈性力學(xué)邊界元法中的近奇異積分可以歸結(jié)為以下兩種形式[32-33]

其中α>0,f(ξ)為規(guī)則函數(shù)(包含雅可比函數(shù)、形函數(shù)等).若采用三點(diǎn)二次單元近似幾何邊界,計(jì)算點(diǎn)P到積分單元的距離函數(shù)r2(ξ)可表示為[32]

其中g(shù)(ξ)0 為規(guī)則函數(shù),a代表計(jì)算點(diǎn)P在積分單元投影的無因次坐標(biāo),b代表計(jì)算點(diǎn)P到積分單元的最短距離.將式(26)代入式(24)和式(25)可得

上述積分在a點(diǎn)被分為兩部分來計(jì)算,可提高數(shù)值計(jì)算精度.在積分區(qū)間[?1,a]和[a,1]上分別引入變量替換x=a?ξ 和x=ξ ?a,則上述近奇異積分可表示為以下形式

其中A=1+a或A=1 ?a.顯然,當(dāng)計(jì)算點(diǎn)P到積分單元的最短距離b很小時,上述積分核函數(shù)在0 點(diǎn)附近會產(chǎn)生劇烈的振蕩特性.為此,引入如下非線性變量替換

其中k=0.5 ln(1+A/b).將變量替換(31)代入式(30)可得

式(29) 的變換類似可得.因g(·)0 為規(guī)則函數(shù),式(32)核函數(shù)分母部分(ek(1+t)?1)2g(t)+11 恒成立,原近奇異性得以消除,可通過傳統(tǒng)高斯積分進(jìn)行精確求解.論文首次嘗試將上述計(jì)算近奇異積分的正則化算法與“特殊裂尖單元”法相結(jié)合,數(shù)值模擬大尺寸比(超薄) 雙材料的界面裂紋問題,下面給出數(shù)值算例.

5 數(shù)值算例

首先,為驗(yàn)證所提“特殊裂尖單元法”求解雙材料界面裂紋問題的有效性,算例1 給出了無限大雙材料板受均勻拉力時斷裂力學(xué)參量和近場物理量的計(jì)算結(jié)果(不涉及近奇異積分的計(jì)算問題).算例2 ～5側(cè)重于大尺寸比(超薄) 雙材料的界面裂紋問題.其中,算例2 和算例3 假設(shè)兩種材料具有相同的彈性模量和泊松比(應(yīng)力強(qiáng)度因子存在精確解).算例4和算例5 討論了不同材料的界面裂紋問題.

例1含有中心裂紋的無限大雙材料板,界面裂紋的長度為2a,裂紋表面受均勻拉力σ0(見圖4).材料1 和材料2 的泊松比為ν1=ν2=0.2.實(shí)際計(jì)算中,令雙材料板的寬度和高度為50a.邊界共劃分280 個二次單元,其中裂紋表面劃分20 個二次單元,裂尖單元相對于裂紋長度為a/10.斷裂力學(xué)參數(shù)和裂尖漸近位移和應(yīng)力場的精確解見文獻(xiàn)[28].

圖4 無限大雙材料板界面裂紋問題Fig.4 An infinite bimaterial plate with an interface crack

當(dāng)雙材料彈性模量的比值μ1/μ2不斷變化時,表1 和表2 給出了復(fù)應(yīng)力強(qiáng)度因子K1和K2的計(jì)算結(jié)果.結(jié)果表明,即使μ1/μ2=1 000,復(fù)應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果與精確解依然十分的吻合.需要注意,不論μ1/μ2如何變化,裂紋表面僅劃分20 個二次單元,且無需對裂尖單元進(jìn)行局部加密.

表1 應(yīng)力強(qiáng)度因子K1/(σ0)的計(jì)算結(jié)果(例1)Table 1 Normalized SIFs K1/(σ0)(Example 1)

表1 應(yīng)力強(qiáng)度因子K1/(σ0)的計(jì)算結(jié)果(例1)Table 1 Normalized SIFs K1/(σ0)(Example 1)

表2 應(yīng)力強(qiáng)度因子K2/(σ0)的計(jì)算結(jié)果(例1)Table 2 Normalized SIFs K2/(σ0√)(Example 1)

圖5 和圖6 分別給出了μ1/μ2=5 時,裂紋表面張開位移和裂尖應(yīng)力場的計(jì)算結(jié)果.結(jié)果表明所提“特殊裂尖單元法”可精確表征界面裂紋漸進(jìn)位移和應(yīng)力場的振蕩特性,在不增加計(jì)算量的前提下,顯著提高斷裂力學(xué)參數(shù)和近場物理量的數(shù)值求解精度.圖7 給出了裂紋表面剖分單元數(shù)量對計(jì)算結(jié)果的影響(μ1/μ2=5).當(dāng)裂紋表面劃分超過20 個單元時,復(fù)應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果趨于穩(wěn)定.

圖5 裂紋表面張開位移(a)與相對誤差(b)Fig.5 Results of the crack-opening-displacements(upper)and relative errors(lower)

圖6 裂紋尖端應(yīng)力場σy 的計(jì)算結(jié)果Fig.6 Results of the near-tip stresses σy

圖7 裂紋表面單元數(shù)量對計(jì)算結(jié)果的影響Fig.7 Variation of the errors with respect to the number of crack elements

例2含有中心裂紋的無限大薄板,裂紋的長度為2a,薄板的寬度為2L,厚度為2h,如圖8 所示.裂紋上下表面受均勻拉力σ.定義薄板的相對厚度為h/L.當(dāng)a/h→∞時,應(yīng)力強(qiáng)度因子的精確解為

其中ν 為材料的泊松比.邊界共劃分126 個二次單元,其中裂紋表面劃分15 個二次單元.實(shí)際計(jì)算中,令L/a=20 來近似無限大板.需要注意,邊界元法僅需在邊界離散,可以從容處理邊界單元的疏密過渡.在模擬薄體結(jié)構(gòu)問題時,不需要根據(jù)薄體結(jié)構(gòu)厚度的變化對網(wǎng)格進(jìn)行加密,這是比有限元模擬薄體結(jié)構(gòu)問題有優(yōu)勢的地方.當(dāng)薄板的相對厚度h/L在10?2～10?9變化時,表3 給出了無量綱應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果.其中,傳統(tǒng)邊界元法(conventional BEM)指對薄體結(jié)構(gòu)中產(chǎn)生的近奇異積分直接采用傳統(tǒng)高斯積分求解.

圖8 含有中心裂紋的薄板Fig.8 A central craked thin beam constrained between two rigid grips

表3 的數(shù)據(jù)表明,當(dāng)薄板的相對厚度h/L<10?4時,傳統(tǒng)邊界元法計(jì)算的應(yīng)力強(qiáng)度因子已經(jīng)完全失效.相比之下,本文算法在薄板的相對厚度小到10?9時,應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果依然與精確解十分吻合.

表3 應(yīng)力強(qiáng)度因子(無量綱K1)的計(jì)算結(jié)果(例2)Table 3 Results of normalized stress intensity factors K1 (Example 2)

例3邊界受剪應(yīng)力的含有中心裂紋的薄板,裂紋的長度為2a,薄板的高度為2H,寬度為2b,如圖9 所示.裂紋長度與板的寬度具有關(guān)系b=4a,薄板的高度設(shè)置為H=15.定義薄板的相對厚度為b/H.邊界共劃分188 個二次單元,其中裂紋表面劃分10 個二次單元.應(yīng)力強(qiáng)度因子的精確解為

圖9 含有中心裂紋的薄板受均勻剪應(yīng)力Fig.9 A central craked thin beam under shear loading

當(dāng)薄板的相對厚度b/H在10?1～10?9變化時,表4 給出了無量綱應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果.表4 的數(shù)據(jù)表明,當(dāng)薄板的相對厚度小于10?4時,傳統(tǒng)邊界元法計(jì)算的應(yīng)力強(qiáng)度因子已經(jīng)完全失效.相比之下,本文算法在薄板的相對厚度小到10?9時,應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果依然與精確解十分吻合.

表4 應(yīng)力強(qiáng)度因子(無量綱K2)的計(jì)算結(jié)果(例3)Table 4 Results of normalized stress intensity factors K2 (Example 3)

例4含有邊裂紋的雙材料板,如圖10 所示.裂紋的長度為a,材料的厚度為h,寬度為L.板的下邊界固定(位移分量為0),左右邊界受拉應(yīng)力,上邊界和裂紋表面不受力.材料1 和材料2 的泊松比為ν1=ν2=0.3,彈性模量μ1/μ21.邊界共劃分320 個二次單元,其中裂紋表面劃分5 個二次單元.定義薄板的相對厚度為h/L,裂紋的相對長度為a/h.固定薄板的厚度為h/L=10?4,裂紋的相對長度在0.01a/h0.99 之間變化.表5 給出了不同材料組合時(彈性模量μ1/μ21),復(fù)應(yīng)力強(qiáng)度因子和的計(jì)算結(jié)果.

圖10 含有邊裂紋的雙材料板(拉應(yīng)力)Fig.10 A bimaterial beam with an edge-crack under uniform tensile stress loading

表5 復(fù)應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果(例4)Table 5 Results of normalized complex stress intensity factors K1/σand K2/σ(Example 4)

表5 復(fù)應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果(例4)Table 5 Results of normalized complex stress intensity factors K1/σand K2/σ(Example 4)

例5如圖11 所示,例5 與例4 具有相同的幾何形狀,只是邊界條件不同.此時,板的下邊界固定(位移分量為零),左右邊界和上邊界不受力,裂紋表面受剪應(yīng)力.同樣,材料1 和材料2 的泊松比為ν1=ν2=0.3,彈性模量μ1/μ21.邊界共劃分320 個二次單元,其中裂紋表面劃分5 個二次單元.與例4類似,固定薄板的厚度為h/L=10?4,裂紋的相對長度在0.01a/h0.99 之間變化.圖12 和圖13 給出了不同材料組合時,復(fù)應(yīng)力強(qiáng)度因子和的計(jì)算結(jié)果.隨著材料參數(shù)μ1/μ2的增大,復(fù)應(yīng)力強(qiáng)度因子的實(shí)數(shù)部分K1不斷增大,而復(fù)應(yīng)力強(qiáng)度因子的復(fù)數(shù)部分K2不斷減小.

圖11 含有邊裂紋的雙材料板(剪應(yīng)力)Fig.11 A bimaterial beam with an edge-crack under uniform shear stress loading

圖12 復(fù)應(yīng)力強(qiáng)度因子K1/(τ)的計(jì)算結(jié)果Fig.12 Results of normalized complex stress intensity factors K1/(τ)

圖13 復(fù)應(yīng)力強(qiáng)度因子K2/(τ)的計(jì)算結(jié)果Fig.13 Results of normalized complex stress intensity factors K2/(τ)

6 結(jié)論

針對求解雙材料界面裂紋復(fù)應(yīng)力強(qiáng)度因子的困難,引入了一種含有復(fù)振蕩因子的新型“特殊裂尖單元”.該方法可直接通過裂尖單元的張開位移或界面應(yīng)力實(shí)現(xiàn)復(fù)應(yīng)力強(qiáng)度因子的精確求解,無需裂尖區(qū)域高密度的網(wǎng)格剖分,也無需借助J 積分、相互作用積分或虛擬裂紋閉合技術(shù)等方法對復(fù)應(yīng)力強(qiáng)度因子進(jìn)行間接求解.該方法也可直接與有限元或其它數(shù)值方法相結(jié)合,用于模擬復(fù)合材料界面裂紋等相關(guān)問題.

此外,將求解薄體結(jié)構(gòu)問題的正則化邊界元法與上述“特殊裂尖單元”法相結(jié)合,成功求解了超薄(大尺寸比) 雙材料的界面裂紋問題.數(shù)值算例表明,所提算法穩(wěn)定,效率高,在不增加計(jì)算量的前提下,顯著提高了裂尖近場力學(xué)參量和斷裂力學(xué)參數(shù)的求解精度和數(shù)值穩(wěn)定性.