蘇海東，祁勇峰，龔亞琦

(長江科學(xué)院材料與結(jié)構(gòu)研究所，武漢 430010)

線彈性斷裂力學(xué)(Linear Elastic Fracture Mechanics，簡稱 LEFM)［1］是斷裂力學(xué)中最早、也是發(fā)展最完善的一個分支。它以線彈性力學(xué)為基礎(chǔ)，將結(jié)構(gòu)視為帶有裂紋的彈性體來研究裂紋的擴(kuò)展問題，其中，裂紋按受力和斷裂特征分為3類:張開型(I型)、滑開型(II型)和撕開型(III型)。根據(jù)經(jīng)典的彈性理論，在裂紋尖端附近的應(yīng)力具有奇異性。Irwin在1957年提出了新的物理量——應(yīng)力強(qiáng)度因子(Stress Intensity Factor，簡稱SIF)來表征裂紋特性，相應(yīng)的有3類應(yīng)力強(qiáng)度因子KⅠ，KⅡ和KⅢ。

計算應(yīng)力強(qiáng)度因子的方法有解析法和數(shù)值方法，前者只適用于簡單形狀。目前基于網(wǎng)格的數(shù)值計算方法主要以有限元法為主，由于奇異性僅存在于裂紋尖端附近很小的局部區(qū)域，往往存在網(wǎng)格劃分困難、收斂慢等問題。近年來出現(xiàn)的數(shù)值流形方法［2］和擴(kuò)展有限元法［3］提出了裂紋在網(wǎng)格內(nèi)穿過的技術(shù)，從而大大降低了網(wǎng)格劃分的難度，后者還引入了裂紋尖端位移場的特殊函數(shù)以加快收斂性，但這些方法仍需進(jìn)一步改進(jìn)。本文基于數(shù)值流形方法，提出裂紋尖端的解析解與其周邊數(shù)值解聯(lián)合應(yīng)用以求解應(yīng)力強(qiáng)度因子的新方法。以下限于平面問題，僅討論Ⅰ型和Ⅱ型裂紋。

1 應(yīng)力強(qiáng)度因子的數(shù)值計算方法及其研究現(xiàn)狀

建立如圖1所示的裂紋尖端處的(r，θ)極坐標(biāo)。以Ⅰ型(張開型)裂紋為例，根據(jù)經(jīng)典的彈性理論，裂紋尖端的位移、應(yīng)力漸近場為:

式中:

計算應(yīng)力強(qiáng)度因子的方法主要有解析法和數(shù)值方法。復(fù)變函數(shù)、積分變換等解析方法一般只能求解簡單形狀的問題，而數(shù)值計算方法適用于復(fù)雜結(jié)構(gòu)，包括目前應(yīng)用最多的有限元法以及新興的數(shù)值流形方法、擴(kuò)展有限元法等。

圖1 裂紋尖端附近應(yīng)力分布Fig.1 Stress distribution near the crack tip

常規(guī)的有限元法要求在網(wǎng)格內(nèi)使用連續(xù)函數(shù)作為插值函數(shù)，且材料性能一致，因此在處理裂紋這樣的強(qiáng)不連續(xù)性問題時，一般要將裂紋尖點(diǎn)設(shè)置為單元結(jié)點(diǎn)，裂紋面設(shè)置為單元的邊，為模擬裂紋尖端的奇異性需要布置非常細(xì)密的網(wǎng)格，而一旦裂紋發(fā)生擴(kuò)展，需沿著裂紋擴(kuò)展方向不斷地修改單元網(wǎng)格和增加新結(jié)點(diǎn)，使用起來很不方便。

數(shù)值流形方法［2］(簡稱流形法)是石根華博士提出的一種新的計算方法，它引入兩套覆蓋(網(wǎng)格)體系:一是用于構(gòu)造物理場近似解的數(shù)學(xué)網(wǎng)格;另一套是表示材料邊界、用于定義積分區(qū)域的物理網(wǎng)格。兩套覆蓋體系相互獨(dú)立，只要求前者在空間上完全包容后者。這樣，材料的物理邊界可以與數(shù)學(xué)網(wǎng)格不匹配，換言之，數(shù)學(xué)網(wǎng)格內(nèi)可以有任意形狀的材料體(稱為流形元)。目前所見的流形法多采用有限元網(wǎng)格作為數(shù)學(xué)網(wǎng)格，在結(jié)點(diǎn)上采用級數(shù)表達(dá)式的覆蓋函數(shù)(一般為多項式級數(shù)，也可以是其他類型級數(shù)甚至是解析解級數(shù)，未知量為級數(shù)的系數(shù))，通過有限元網(wǎng)格的形函數(shù)(流形法中稱為權(quán)函數(shù))相連形成網(wǎng)格內(nèi)的插值函數(shù)。當(dāng)采用多項式覆蓋函數(shù)時，一般應(yīng)用單純形精確積分法以保證數(shù)學(xué)網(wǎng)格內(nèi)任意形狀流形元的積分精度。為避免歧義，以下所述的“網(wǎng)格”均指數(shù)學(xué)網(wǎng)格。

對于不連續(xù)分析問題，流形法引入不連續(xù)覆蓋來模擬裂紋［2］，如圖2所示的具有1條裂紋的材料體，裂紋將網(wǎng)格分割成不連續(xù)的區(qū)域，每個區(qū)域稱為一個物理覆蓋，在裂紋的兩邊采用不同的覆蓋編號，比如，結(jié)點(diǎn)覆蓋7被裂紋分割成71和72，在網(wǎng)格7-8-13-12中，材料體被裂紋分割成兩個流形元，裂紋下邊的流形元覆蓋為71-81-131-121，裂紋上邊為72-82-132-122。這樣，裂紋可以在網(wǎng)格內(nèi)部穿過，巧妙地解決了常規(guī)有限元法裂紋面必須與單元邊一致、裂紋擴(kuò)展后需要重新劃分網(wǎng)格的問題。但目前，流形法還難以處理裂紋尖端停留在網(wǎng)格中間的情況，如圖2中的網(wǎng)格8-9-14-13。

圖2 有一條裂紋的材料體的矩形網(wǎng)格覆蓋Fig.2 Rectangular meshes covering material domain with one crack

擴(kuò)展有限元法［3］借鑒了流形法思想，裂紋同樣可以在有限單元中穿過，對于裂紋尖端停留在網(wǎng)格中間的情況，在網(wǎng)格內(nèi)引入附加基函數(shù)來擴(kuò)展自由度，基函數(shù)一般取為裂紋尖端漸近位移場(見式(1))中的等以反映裂紋尖端的奇異性。反過來，也有流形法的相關(guān)文獻(xiàn)［4］借鑒擴(kuò)展有限元的這種方式處理裂紋尖端停留在網(wǎng)格中間的情況。

以上方法求解應(yīng)力強(qiáng)度因子時，都需要通過計算裂紋尖端附近的位移或應(yīng)力來推求應(yīng)力強(qiáng)度因子，對位移或應(yīng)力的精度要求高。實(shí)際上，有一種更直接的方法，利用裂紋尖端附近的Williams位移解析級數(shù)［5］直接求得應(yīng)力強(qiáng)度因子。由于引入了與位移自由度不同的級數(shù)系數(shù)作為未知數(shù)，因此裂紋尖端附近的這種單元被稱為雜交裂縫單元(Hybrid Crack Element，簡稱HCE)。該單元與周邊以位移為自由度的常規(guī)單元的連接問題是該方法的研究重點(diǎn)，在常規(guī)有限元法的框架下，這種連接過渡是比較困難的，一般要采用特殊的網(wǎng)格劃分［6］，或者如文獻(xiàn)［7］提出的廣義參數(shù)單元。

流形法具有解析解和數(shù)值解聯(lián)合應(yīng)用的方便性，作者曾在文獻(xiàn)［8］中應(yīng)用這種特性很好地求解了帶有特殊無限域流場的流固耦合問題。本文沿用這種思路，運(yùn)用裂紋尖端附近的Williams解析解與周邊數(shù)值解聯(lián)合求解應(yīng)力強(qiáng)度因子。

2 裂紋尖端位移場的Williams級數(shù)解析解及其覆蓋函數(shù)

如圖1所示的裂紋尖端區(qū)域，Williams位移級數(shù)解為(裂紋邊界滿足無應(yīng)力條件)

式中:

ai，bi為待求系數(shù)，其中a0，b0代表剛體位移項。應(yīng)力強(qiáng)度因子與ai，bi的關(guān)系是

顯而易見，Ⅰ型裂紋表達(dá)式(1)就是式(3)中取i=1，m=1及 bi=0的特殊情況。

如圖3所示，在裂紋尖端所在的矩形網(wǎng)格6-7-11-10內(nèi)，每個結(jié)點(diǎn)均采用式(3)的覆蓋函數(shù)形式，并用矩形網(wǎng)格的形函數(shù)wj連接，形成網(wǎng)格內(nèi)的插值函數(shù)，寫成矩陣形式:

本文提出了采用結(jié)點(diǎn)自由度之間的強(qiáng)制約束方法，將該網(wǎng)格的結(jié)點(diǎn)7，11，10約束到結(jié)點(diǎn)6，即令，則根據(jù)形函數(shù)的關(guān)系，從而在該網(wǎng)格內(nèi)形成獨(dú)立的解析解覆蓋

正好就是式(3)的矩陣形式。

以下推導(dǎo)應(yīng)變矩陣和剛度矩陣。根據(jù)式(5)，應(yīng)變子矩陣的第i項為

其中:

圖3 裂紋尖端附近的矩形網(wǎng)格Fig.3 Rectangular meshes around the crack tip

另外:(x0，y0)為裂紋尖端坐標(biāo)。

剛度子矩陣

式中:［D］為彈性矩陣，A為流形元面積。

3 裂紋周邊網(wǎng)格內(nèi)的解析解與數(shù)值解的覆蓋聯(lián)系

除了裂紋尖端所在的網(wǎng)格外，在裂紋周邊與該網(wǎng)格相連的其他網(wǎng)格內(nèi)，位移統(tǒng)一表示為

其中:cnk，dnk為待求的多項式系數(shù);n，k為多項式階次(h為多項式的最高階次);wj，wl為相應(yīng)于各結(jié)點(diǎn)的形函數(shù);j和l的個數(shù)根據(jù)網(wǎng)格而定，但保持j+l=4。顯然，取l=0則退化到解析解覆蓋式(6)。

如圖3所示裂紋周邊網(wǎng)格，比如:網(wǎng)格3-4-8-7中的結(jié)點(diǎn)7為解析解覆蓋，其它3個結(jié)點(diǎn)為數(shù)值解覆蓋;網(wǎng)格7-8-12-11中的結(jié)點(diǎn)7和結(jié)點(diǎn)11為解析解覆蓋(強(qiáng)制約束為同一覆蓋函數(shù))，結(jié)點(diǎn)8和結(jié)點(diǎn)12為數(shù)值解覆蓋;網(wǎng)格5-6-10-9采用不連續(xù)的覆蓋，實(shí)現(xiàn)裂紋兩邊的獨(dú)立運(yùn)動。

對于解析解覆蓋結(jié)點(diǎn)j處的第i項，其應(yīng)變矩陣為

而對于數(shù)值解覆蓋結(jié)點(diǎn)l處，第q項的應(yīng)變矩陣為

解析覆蓋與數(shù)值覆蓋相關(guān)的剛度子矩陣為

由式(7)至式(13)可見，剛度矩陣積分中具有非多項式的函數(shù)，因此基于多形式函數(shù)的單純形積分公式無法應(yīng)用。本文采用三角形區(qū)域的Hammer積分［9］，為保障數(shù)值積分的精度，考慮一種逐步的積分區(qū)域加密方式，如圖2中的網(wǎng)格8-9-14-13中形成的流形元(見圖4)，取其中心與流形元各邊相連形成三角形積分區(qū)域即可進(jìn)行Hammer數(shù)值積分，若積分精度不夠，則在各三角形中，將邊中點(diǎn)相連分解成4個小三角形，以此類推。一般而言，只需加密1次、最多2次就可以滿足積分精度要求。

圖4 積分區(qū)域加密Fig.4 Subdivided integration domains

4 算例

以無限長條為例，考慮以下幾種典型裂紋。

4.1 兩端受均布拉力的內(nèi)部裂紋

如圖5(a)所示，在寬度為B的無限長條中央有一條長度為2a的內(nèi)部裂紋，長條兩端受均布拉力p的作用。當(dāng)4a＜B時，應(yīng)力強(qiáng)度因子理論值［10］為取 a=0.05 m，B=0.8 m，p=0.3 kN/cm，則 KⅠ理論值為1.20。

4.2 兩端受均布拉力的邊界裂紋

見圖5(b)，長條的寬度為w=0.4 m，其他同圖5(a)，應(yīng)力強(qiáng)度因子理論值［10］為其中，f(a/w)=1.12-0.23(a/w)+10.6(a/w)2-21.7(a/w)3+30.4(a/w)4，該算例中 KⅠ數(shù)值為1.45。

4.3 在裂紋嘴處受一對集中切向力作用的邊界裂紋

如圖5(c)所示，無限長條在其裂紋嘴處有一對集中切向力F=3 kN，其應(yīng)力強(qiáng)度因子理論值［10］為，其中

該算例中KⅡ理論值為1.98。

流形元網(wǎng)格見圖6，圖6(a)為整體網(wǎng)格，圖6(b)為裂紋附近的網(wǎng)格，網(wǎng)格尺寸為0.012 5～0.02 m，其中對于本文4.1節(jié)裂紋，考慮對稱性取一半計算，對稱面施加法向約束，圖中橫向粗線表示裂紋，豎向粗線為對稱面，可見流形網(wǎng)格與物理邊界不一致。而對于本文4.2節(jié)和4.3節(jié)的裂紋，在側(cè)向不加法向約束。

圖6 流形元網(wǎng)格Fig.6 NMM meshes

對于兩端受均布拉力的內(nèi)部裂紋(本文4.1節(jié))，應(yīng)力強(qiáng)度因子KⅠ的計算結(jié)果見表1，隨著Williams級數(shù)的階數(shù)m以及周邊數(shù)值解多項式函數(shù)最高階數(shù) n的增加，KⅠ的漸近值為1.23，與理論值1.20很接近。

表1 應(yīng)力強(qiáng)度因子KⅠ(兩端受均布拉力的內(nèi)部裂紋情況)Table 1 Stress Intensity Factor KⅠ(Two ends subjected to uniform tension，internal crack)

計算兩端受均布拉力的邊界裂紋(本文4.2節(jié))，理論值為1.45，計算得到的 KⅠ漸近值也為1.45。

再計算邊界裂紋受一對集中切向力的裂紋(本文4.3節(jié))，理論值為1.98，計算得到的 KⅡ漸近值為2.00。

從以上計算結(jié)果可見，本文方法可以得到精度很高的應(yīng)力強(qiáng)度因子計算值。以下對裂紋尖端附近的網(wǎng)格密度問題進(jìn)行討論。

如圖7所示，將數(shù)學(xué)網(wǎng)格的尺寸加大，裂紋尖端附近的網(wǎng)格尺寸為0.03 ～0.05 m。對于兩端受均布拉力的內(nèi)部裂紋，應(yīng)力強(qiáng)度因子KⅠ的計算結(jié)果見表2，可見計算精度明顯降低，m至少取到7階以上，KⅠ才能達(dá)到1左右。隨著階數(shù)的升高，計算也表現(xiàn)出一定的不穩(wěn)定性，這表明本文方法對裂紋附近的網(wǎng)格密度仍有一定的要求，當(dāng)網(wǎng)格密度不夠時，計算結(jié)果不理想。至于如何判斷是否達(dá)到理想的網(wǎng)格密度，這涉及到流形法的精度分析問題，需要更深入的研究。

圖7 裂紋周邊的局部流形元網(wǎng)格(粗網(wǎng)格)Fig.7 Coarse NMM meshes around the crack tip

表2 應(yīng)力強(qiáng)度因子KⅠ(兩端受均布拉力的內(nèi)部裂紋，粗網(wǎng)格)Table 2 Stress intensity factor KⅠ(Two ends subjected to uniform tension，internal crack，coarse meshes)

5 結(jié)語

由于裂紋尖端位移和應(yīng)力分布的復(fù)雜性，用常規(guī)有限元網(wǎng)格的插值函數(shù)去逼近裂紋尖端的位移和應(yīng)力是不易做到的，必須在裂紋尖端附近布置非常細(xì)密的網(wǎng)格。近年來發(fā)展起來的擴(kuò)展有限元采用了裂紋尖端漸近位移場中的特征函數(shù)，能夠比常規(guī)有限元插值函數(shù)更好地反映出裂紋尖端的奇異性，因而收斂較快。但擴(kuò)展有限元只是使用了裂紋尖端漸近位移場的部分特征函數(shù)，而本文方法采用的Williams解析級數(shù)是對裂紋尖端位移場的最佳逼近，同時，周邊數(shù)值解采用的高階多項式插值函數(shù)相對于常規(guī)有限元和一般擴(kuò)展有限元的普通插值多項式而言逼近的精度更高，因此本文方法收斂最快。

常規(guī)有限元和擴(kuò)展有限元計算應(yīng)力強(qiáng)度因子通常采用的“直接”法，需要通過裂紋尖端附近的位移或應(yīng)力去擬合裂紋尖端的表達(dá)式進(jìn)而推求應(yīng)力強(qiáng)度因子;而所謂J積分的“間接”法，也需要在裂紋上、下表面定義的一個回線上進(jìn)行積分，不同的回線得到的結(jié)果也會有差異。這些操作會引入額外誤差。本文方法將應(yīng)力強(qiáng)度因子作為方程組的未知數(shù)，不需要其他操作，是最“直接”最便利的方法，不會引入額外誤差。

對于雜交裂縫單元(HCE)方法，雖然也是利用裂紋尖端附近的Williams解析解級數(shù)直接求應(yīng)力強(qiáng)度因子，但這種單元與周邊以位移為自由度的常規(guī)單元的連接過渡是比較困難的。在這一點(diǎn)上，本文方法具有較大的優(yōu)勢，流形法的不連續(xù)覆蓋技術(shù)可以使裂紋在網(wǎng)格內(nèi)穿過，這是常規(guī)有限元法以及在其框架下的HCE方法很難做到的。

綜上所述，本文方法相比現(xiàn)有方法具有更快速的收斂性和更大的方便性。另外，對于流形法的研究而言，本文方法的意義在于:

(1)雖然流形法自誕生之時就以不連續(xù)覆蓋的非連續(xù)變形分析方法而聞名于世，但如何處理裂紋尖端停留在網(wǎng)格內(nèi)部的情況一直是其難以解決的問題，甚至落后于較晚出現(xiàn)的擴(kuò)展有限元法，而本文方法成功解決了此問題。

(2)本文方法是繼文獻(xiàn)［8］之后的數(shù)值解和解析解聯(lián)合運(yùn)用的又一典型實(shí)例，再一次驗證了梁國平教授在文獻(xiàn)［2］的序中所述的“采用數(shù)值流形方法很容易就解決了人們多年來研究的至今尚未有好的解決辦法的局部區(qū)域解析法以及有限元與解析法相結(jié)合的方法”，充分體現(xiàn)出流形法的優(yōu)勢。

本文方法雖然有更快的收斂性，但對裂紋尖端附近的網(wǎng)格密度仍有一定的要求，必要時需要加密網(wǎng)格。采用常規(guī)的網(wǎng)格加密方法，從普通大小的網(wǎng)格過渡到裂紋尖端處小一至兩個數(shù)量級尺寸的網(wǎng)格是比較困難的。我們最近研究的部分重疊覆蓋的流形法，可以方便地做到大小網(wǎng)格之間的協(xié)調(diào)過渡，非常適合于裂紋尖端的網(wǎng)格加密，將在另文介紹。

致謝:感謝數(shù)值流形方法的發(fā)明者石根華博士的指導(dǎo)。

［1］范天佑.斷裂理論基礎(chǔ)［M］.北京:科學(xué)出版社，2003.(FAN Tian-you.Basic Principle of Fracture［M］.Beijing:Science Press，2003.(in Chinese))

［2］石根華.數(shù)值流形方法與非連續(xù)變形分析［M］.裴覺民譯.北京:清華大學(xué)出版社，1997.(SHI Gen-hua.Numerical Manifold Method(NMM)and Discontinuous Deformation Analysis(DDA)［M］.Translated by PEI Jue-min.Beijing:Tsinghua University Press，1997.(in Chinese))

［3］李錄賢，王鐵軍.擴(kuò)展有限元法(XFEM)及其應(yīng)用［J］.力學(xué)進(jìn)展，2005，35(1):5-20.(LI Lu-xian，WANG Tie-jun.The Extended Finite Element Method and Its Applications-A Review［J］.Advances in Mechanics，2005，35(1):5-20.(in Chinese))

［4］李樹忱，程玉民.考慮裂紋尖端場的數(shù)值流形方法［J］.土木工程學(xué)報，2005，38(7):96-101.(LI Shuchen，CHENG Yu-min.Numerical Manifold Method for Crack Tip Fields［J］.China Civil Engineering Journal，2005，38(7):96-101.(in Chinese))

［5］WILLIAMS M L.On the Stress Distribution at the Base of a Stationary Crack［J］.Journal of Applied Mechanics，1957，24:109-114.

［6］KARIHALOO B L，XIAO Q Z.Implementation of HCE on a General FE Mesh for Interacting Multiple Cracks［C］∥Proceedings of the European Congress on Computational Methods in Applied Sciences and Engineering，ECCOMAS，Jyv?skyl?:July 24-28，2004.

［7］彭 俚.結(jié)構(gòu)斷裂分析的Williams廣義參數(shù)單元［D］.南寧:廣西大學(xué)，2008.(PENG Li.Williams Element with Generalized DOF for Structural Fracture Analysis［D］.Nanning:Guangxi University，2008.(in Chinese))

［8］蘇海東，黃玉盈.數(shù)值流形方法在流固耦合諧振分析中的應(yīng)用［J］.計算力學(xué)學(xué)報，2007，24(6):823-828.(SU Hai-dong，HUANG Yu-ying.Application of Numerical Manifold Method in Fluid-Solid Interaction Harmonic Analysis［J］.Chinese Journal of Computational Mechanics，2007，24(6):823-828.(in Chinese))

［9］王勖成，邵 敏.有限單元法的基本概念和數(shù)值方法［M］.北京:清華大學(xué)出版社，1988.(WANG Xu-cheng，SHAO Min.Basic Principle of the FEM and Numerical Method［M］.Beijing:Tsinghua University Press，1988.(in Chinese))

［10］中國航空研究院.應(yīng)力強(qiáng)度因子手冊［M］.北京:科學(xué)出版社，1993.(Chinese Institute of Aviation Research.Manual of Stress Intensity Factors［M］.Beijing:Science Press，1993.(in Chinese))