劉丹丹，李曉川

(沈陽(yáng)工業(yè)大學(xué)建筑工程學(xué)院，遼寧 沈陽(yáng) 110870)

從20世紀(jì)50年代開始，固體力學(xué)的一個(gè)分支斷裂力學(xué)快速發(fā)展了起來(lái)，并為現(xiàn)代化的生產(chǎn)提供了理論依據(jù)。半個(gè)多世紀(jì)過去了，斷裂力學(xué)[1]這一領(lǐng)域的知識(shí)理論已經(jīng)發(fā)展的相當(dāng)成熟。斷裂力學(xué)有助于人們認(rèn)識(shí)和理解材料中裂紋的產(chǎn)生、擴(kuò)展和最終破壞的過程，并為結(jié)構(gòu)損傷容限設(shè)計(jì)提供理論基礎(chǔ)。

斷裂力學(xué)中的3個(gè)基本斷裂參數(shù)是應(yīng)力強(qiáng)度因子、路徑無(wú)關(guān)積分（J積分）和應(yīng)變能釋放率。應(yīng)力強(qiáng)度因子的計(jì)算主要依賴于裂紋前端的局部應(yīng)力場(chǎng)，它描述了彈性裂紋尖端應(yīng)力場(chǎng)的強(qiáng)弱。J積分描述的是由于裂紋的存在所吸收的能量，而能量釋放率是描述產(chǎn)生新的裂紋面所需要的能量。這兩個(gè)基于能量的參數(shù)是等價(jià)[2]的。隨著計(jì)算機(jī)硬件和軟件的迅速發(fā)展，用數(shù)值方法計(jì)算斷裂參數(shù)就變的切實(shí)可行。很多數(shù)值方法被嘗試應(yīng)用于斷裂參數(shù)的計(jì)算，如有限差分法、邊界元法和發(fā)展起來(lái)的無(wú)網(wǎng)格法，但是由于缺少軟件的支持，這些方法應(yīng)用實(shí)例相對(duì)缺乏。而有限元法成功的應(yīng)用于很多工業(yè)部門，并發(fā)展了針對(duì)三個(gè)基本斷裂參數(shù)的數(shù)值計(jì)算技術(shù)。用應(yīng)力或應(yīng)變外推法計(jì)算應(yīng)力強(qiáng)度因子，等效積分區(qū)域技術(shù)計(jì)算J積分，全局或者局部的虛擬裂紋擴(kuò)展技術(shù)和虛擬裂紋閉合法計(jì)算應(yīng)變能釋放率。

1 外推法

1.1 基于應(yīng)力的外推法

對(duì)于線彈性材料，應(yīng)力在裂紋尖端處是趨于無(wú)窮大的，也就是應(yīng)力奇異性。所以網(wǎng)格尺寸越小，裂紋前端單元內(nèi)的積分點(diǎn)的應(yīng)力值就會(huì)越大。從這可以知道，應(yīng)力值對(duì)于網(wǎng)格的尺寸是不收斂的。于是引入了應(yīng)力強(qiáng)度因子的概念來(lái)克服由于應(yīng)力奇異帶來(lái)的數(shù)學(xué)上的困擾，用來(lái)描述裂紋尖端附近應(yīng)力的強(qiáng)弱。外推法是借助于有限元的方法來(lái)計(jì)算應(yīng)力強(qiáng)度因子，在有限元分析中，裂尖前面的單元里積分點(diǎn)上的應(yīng)力值以及其對(duì)應(yīng)的積分點(diǎn)坐標(biāo)值是很容易就可以直接讀出的。雖然無(wú)法用數(shù)值方法直接計(jì)算裂尖處的應(yīng)力強(qiáng)度因子，可是裂尖前端那些非奇異的應(yīng)力值卻是可知的。對(duì)于每一個(gè)ri（距離裂紋尖端的極半徑）大于零都有一個(gè)非奇的應(yīng)力值σyi和對(duì)應(yīng)的KIi：KIi=σyi，然后可以構(gòu)造數(shù)據(jù)對(duì)（ri,KI）i，用最小二乘法擬合數(shù)據(jù)點(diǎn)。在這里假定了兩者之間線性關(guān)系為：

上式中的N 是外推法所采用的數(shù)據(jù)對(duì)象，其截距的物理意義就是所要計(jì)算的應(yīng)力強(qiáng)度因子。

可是這個(gè)方法存在幾個(gè)問題。第一，沒有明確規(guī)定采取多少對(duì)數(shù)據(jù)點(diǎn)進(jìn)行外推計(jì)算得到最后的精度高。一般情況下是采用許多點(diǎn)直至KI對(duì)N 收斂于某個(gè)值才行。第二，靠近裂尖處的數(shù)據(jù)反而是引起誤差的主要來(lái)源，而相反的距離裂尖越遠(yuǎn)越滿足線性的假設(shè)，這樣使得外推法所采用的點(diǎn)往往不包括貼近裂尖的點(diǎn)，說(shuō)明線性外推對(duì)裂尖來(lái)說(shuō)并不很精確，而在斷裂分析中又恰好最應(yīng)該多關(guān)注裂紋尖端，因?yàn)閼?yīng)力強(qiáng)度因子的作用就是描述裂尖應(yīng)力信息，計(jì)算一個(gè)用來(lái)描述裂尖應(yīng)力場(chǎng)的量所采用的方法卻不能依靠裂尖信息，而必須用遠(yuǎn)場(chǎng)數(shù)據(jù)來(lái)計(jì)算，這在邏輯上存在悖論。

1.2 基于應(yīng)變的外推法

在很多的商業(yè)軟件中，位移是求解的基本變量，而應(yīng)力是通過應(yīng)變和位移聯(lián)系起來(lái)是次要的變量，所以位移的外推法的計(jì)算精度要高。對(duì)于確定的距離裂紋尖端r 處，裂紋后端垂直位移v 的數(shù)據(jù)可以在有限元分析中直接讀取。和基于應(yīng)力的外推法一樣，也可以構(gòu)造數(shù)據(jù)對(duì)（ri,KIi）：

然后利用最小二乘法來(lái)擬合數(shù)據(jù)點(diǎn)，得到相應(yīng)的應(yīng)力強(qiáng)度因子數(shù)值。同樣的，基于應(yīng)力的外推法存在的問題也存在于此方法中。如果想解決問題，使用多項(xiàng)式外推應(yīng)該是最直接的方法，可是并不是所有的情況都存在確定的多項(xiàng)式。因此，引入奇異單元和疊合單元[3]理論上在裂尖處采用這些單元，在其他區(qū)域采用常規(guī)的單元，可以極大地減少網(wǎng)格的數(shù)量同時(shí)保持計(jì)算精度。然而采用這些單元在使用時(shí)非常的不方便，尤其是對(duì)于裂紋擴(kuò)展的問題，而且適用范圍也很有限，無(wú)法滿足實(shí)際情況多樣性和復(fù)雜性的需求。

2 等效積分區(qū)域法

J積分是由Rice[4]提出來(lái)的一個(gè)處理非線性斷裂問題的斷裂參數(shù)，這個(gè)參數(shù)的引入是基于能量守恒的概念，因此對(duì)裂紋尖端應(yīng)力奇異性的依賴程度比較弱，不需要對(duì)裂紋尖端處的單元進(jìn)行特殊的處理，而是從全局的角度分析處理斷裂問題。J積分的數(shù)學(xué)表達(dá)式為：

上式中，ui是位移矢量的分量；ds是沿著積分路徑的微小增量；w 是應(yīng)變能密度因子，其定義是：

上式中σij和εij分別是應(yīng)力張量和應(yīng)變張量。

由于Rice證明了J積分的數(shù)值不隨著積分回路的改變而改變，因此它也叫做路徑的無(wú)關(guān)積分。在這里可以看到J積分的表達(dá)式中，只有應(yīng)力、應(yīng)變和位移三場(chǎng)沒有材料屬性。但是在計(jì)算應(yīng)力的時(shí)候是需要用到材料各向同性或者異性的本構(gòu)關(guān)系。所以，這個(gè)表達(dá)式只是使用于一般的熱彈塑性和各向異性材料的問題。實(shí)際上這個(gè)表達(dá)式并不適用于數(shù)值計(jì)算，為此Shih和Raju[5]等人提出了等效積分區(qū)域法來(lái)達(dá)到對(duì)J積分進(jìn)行數(shù)值計(jì)算?；谟昧鸭飧浇囊粋€(gè)有限區(qū)域來(lái)代替積分回路進(jìn)行J積分計(jì)算的思想，表達(dá)式可以轉(zhuǎn)化為：

上式中δij是Kronic函數(shù)，q 是輔助函數(shù)，它只是一個(gè)數(shù)學(xué)處理，目的是為了讓積分表達(dá)式更利于采用數(shù)值計(jì)算方法。J積分是遠(yuǎn)場(chǎng)回路積分，它具有較好的精度，也不需要對(duì)裂尖處的單元進(jìn)行特殊處理，不過由于沒有分離斷裂模式，對(duì)于混合斷裂型問題而言，還需要額外處理才能進(jìn)行完整的分析。就算對(duì)于平面問題，它的計(jì)算過程都顯得相當(dāng)?shù)膹?fù)雜。所以盡管它能夠用于裂紋擴(kuò)展問題但是實(shí)際上操作很不方便，而且在有限元的分析中實(shí)用性很低。

3 虛擬裂紋擴(kuò)展法

3.1 全局的虛擬裂紋擴(kuò)展法

Irwin[6]提出了應(yīng)變能釋放率的概念。考慮有一個(gè)厚度為B的二維裂紋體，其裂紋長(zhǎng)度是a，那么應(yīng)變能釋放率就是產(chǎn)生面積為ΔA 的新裂紋面所需要的能量，于是：

上式中∏=U-W 是裂紋體的勢(shì)能，U 是裂紋體的應(yīng)變能，W是外力功，Δa 是裂紋的微小擴(kuò)展量。

從上式可以看出，應(yīng)變能釋放率的計(jì)算要求Δa 趨近于零，很明顯對(duì)于有限元分析這種數(shù)值方法來(lái)說(shuō)，無(wú)法達(dá)到這個(gè)極限。所以，采用兩步分析過程的虛擬裂紋擴(kuò)展法。在第一步的分析中，通過有限元分析獲得裂紋長(zhǎng)度是a 的裂紋體的勢(shì)能（∏1=U1-W1）；第二步分析中，通過有限元分析獲得裂紋長(zhǎng)度是a+Δa 的裂紋體的勢(shì)能（∏2=U2-W2）。若和網(wǎng)格尺寸相關(guān)的Δa足夠小，那么就可以很好的將應(yīng)變能釋放率近似為：

全域虛擬裂紋擴(kuò)展法有一個(gè)局限性就是只能夠得到總的應(yīng)變能釋放率，不能分離斷裂模式。

3.2 局部的虛擬裂紋擴(kuò)展法

基于將裂紋閉合一個(gè)微小的擴(kuò)展增量所需要的功和勢(shì)能的改變是等效的這一思想，可以用裂紋的閉合積分來(lái)計(jì)算裂尖的能量釋放率，則局部的形式為：

其中應(yīng)變能分量為：

假設(shè)裂紋的方向是沿著X 軸方向，σyy是張開斷裂模式中沿著閉合裂紋面上的法向應(yīng)力，txy是滑移模式中沿著閉合裂紋面上的切向應(yīng)力。Δu 和Δv 是閉合裂紋張開后裂尖后面閉合裂紋面上的位移分量。上兩式可以用兩個(gè)步驟的分析過程計(jì)算出來(lái)。如果網(wǎng)格尺寸相關(guān)的Δa 很小，就可以將應(yīng)變能釋放分量近似為：

如果直接用上兩式計(jì)算應(yīng)變能釋放率分量則需要節(jié)點(diǎn)上的應(yīng)力值，因?yàn)楣街惺茄刂]合裂紋線對(duì)應(yīng)力進(jìn)行數(shù)值積分而閉合裂紋線一般位于有限元單元的邊上。在有限元分析中，當(dāng)把應(yīng)力外推到單元的節(jié)點(diǎn)上，或者當(dāng)單元接近裂紋尖端時(shí)，所得到的應(yīng)力時(shí)非常不精確的。為了避免使用不精確的應(yīng)力值，采用節(jié)點(diǎn)力來(lái)代替應(yīng)力的積分。用節(jié)點(diǎn)力和節(jié)點(diǎn)位移計(jì)算出的應(yīng)變能釋放率為：

用節(jié)點(diǎn)力和節(jié)點(diǎn)位移計(jì)算應(yīng)變能釋放率有很多的優(yōu)點(diǎn)是：①公式中只有節(jié)點(diǎn)力和節(jié)點(diǎn)位移這兩個(gè)基本變量，任何有限元計(jì)算的商業(yè)軟件都可以直接輸出；②公式不需要對(duì)應(yīng)力進(jìn)行數(shù)值積分，這樣減少了計(jì)算量；③有限元網(wǎng)格尺寸的影響小，為了保證近似精度，需要采用合理的網(wǎng)格密度，但是用比較粗糙的網(wǎng)格也能獲得令人滿意的結(jié)果；④無(wú)論是線性的還是非線性的材料都可以采用這個(gè)方法。

而此方法的缺點(diǎn)是分析過程必須是2步，且裂紋長(zhǎng)度在這兩步中是不同的，這樣增加了網(wǎng)格準(zhǔn)備工作量，尤其對(duì)于三維問題時(shí)，同時(shí)，對(duì)于裂紋擴(kuò)展問題，由于需要用上一步的計(jì)算結(jié)果來(lái)不斷的準(zhǔn)備新的網(wǎng)格，使得研究起來(lái)相當(dāng)?shù)牟环奖恪?/p>

4 虛擬裂紋閉合法

為了避免虛擬裂紋擴(kuò)展法的不便之處。Rybicki和Kanninen提出[7]適用于二維平面問題的一步分析法即修正的裂紋閉合積分，也就是后來(lái)的虛擬裂紋閉合法。

虛擬裂紋閉合法的基本假設(shè)是虛擬裂尖和實(shí)際裂尖后面張開位移近似相同。于是，公式轉(zhuǎn)化為：

其中，F(xiàn)y1和Fx1為裂紋尖端節(jié)點(diǎn)Y 和X 方向的節(jié)點(diǎn)力；Δu3,4和Δv3,4為緊挨裂紋尖端后面節(jié)點(diǎn)在X 和Y 方向的相對(duì)位移；B 為裂紋體的厚度。Δa 為裂紋尖端前面的虛擬裂紋擴(kuò)展量。

虛擬裂紋閉合法在保留了虛擬裂紋擴(kuò)展法的所有優(yōu)點(diǎn)的同時(shí)，克服了它的所有缺點(diǎn)，只要求一步的有限元分析。這使得此方法有很大的吸引力，尤其是那些用有限元計(jì)算商業(yè)軟件對(duì)斷裂問題進(jìn)行分析研究的工程師們。

5 結(jié) 論

由上述分析可知，在結(jié)合有限元分析斷裂問題的數(shù)值計(jì)算方法中，虛擬裂紋閉合法非常簡(jiǎn)便高效。

[1]丁遂棟.斷裂力學(xué)[M].北京:機(jī)械工業(yè)出版社,1991.

[2]沈成康.斷裂力學(xué)[M].上海：同濟(jì)大學(xué)出版社,1996.

[3]Blandford,G.E..Two-dimensional stress intensity factor computations using the boundary element method.International Journal for Numerical Methodsin Engineering,1981.

[4]Rice,J.R..APath Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks,Journal of Applied Mechanics,1968.

[5]Moura,B.,and Shih,C.F..A Treatment of Crack Tip Contour Integrals.International Journal of Fracture,1987.

[6]Irwin GR.Crack-extension forcefor a part-through crack in aplate.J Appl Mech,1962.

[7]Rybicki EF,Kanninen MF.A finite element calculation of stress intensity factors by a modified crack closure integral.Engineering Fracture Mechanics,1977(9).