祁勇峰，蘇海東，龔亞琦

(長(zhǎng)江科學(xué)院 材料與結(jié)構(gòu)研究所；水利部水工程安全和病害防治工程技術(shù)研究中心，武漢 430010)

應(yīng)力強(qiáng)度因子準(zhǔn)則[1]是目前最常用的結(jié)構(gòu)裂縫擴(kuò)展準(zhǔn)則之一，基于線彈性斷裂力學(xué)的應(yīng)力強(qiáng)度因子(Stress Intensity Factor，簡(jiǎn)稱SIF)是用來表征裂紋尖端附近應(yīng)力場(chǎng)和應(yīng)變場(chǎng)強(qiáng)度，是控制裂紋尖端應(yīng)力場(chǎng)和應(yīng)變場(chǎng)強(qiáng)度的關(guān)鍵參數(shù)，在裂紋擴(kuò)展分析中有著極其重要的地位。由于應(yīng)力強(qiáng)度因子取決于外力的大小和分布、結(jié)構(gòu)的幾何條件以及裂縫的形狀和位置，實(shí)際上只有少數(shù)問題存在解析解，對(duì)于復(fù)雜幾何形狀和加載條件的問題，只能通過數(shù)值方法來計(jì)算。目前裂縫擴(kuò)展分析的主流數(shù)值計(jì)算方法有有限元法、擴(kuò)展有限元法、數(shù)值流形法等。

有限元法[2-4]是目前分析裂縫等不連續(xù)問題的主要方法，為體現(xiàn)裂紋尖端(下簡(jiǎn)稱裂尖)周圍的應(yīng)力集中和奇異性，往往需要在裂尖附近的復(fù)雜應(yīng)力區(qū)布置高密度的單元網(wǎng)格，導(dǎo)致單元數(shù)目非常龐大；另外，在模擬裂縫擴(kuò)展過程時(shí)，需不斷重構(gòu)網(wǎng)格，因此，有限元法對(duì)網(wǎng)格的要求和依賴性極大地影響了計(jì)算效率。擴(kuò)展有限元法[5-7]通過在單元插值函數(shù)中引入能夠反映裂縫面特性的不連續(xù)階躍函數(shù)及反映裂尖局部特性的裂尖漸進(jìn)位移場(chǎng)函數(shù)，裂縫可以穿過單元內(nèi)部，裂縫擴(kuò)展以后無需重新劃分單元網(wǎng)格，采用同一網(wǎng)格就可以分析任意位置的裂縫問題，克服了常規(guī)有限元進(jìn)行裂縫擴(kuò)展分析的缺點(diǎn)，極大地簡(jiǎn)化了前處理工作。數(shù)值流形方法[8-11]引入不連續(xù)覆蓋模擬裂縫，裂縫可以在網(wǎng)格內(nèi)部穿過，巧妙地解決了常規(guī)有限元法裂縫面必須與單元邊一致、裂縫擴(kuò)展后需要重新劃分網(wǎng)格的問題。相對(duì)于擴(kuò)展有限元法設(shè)立不連續(xù)階躍函數(shù)的方式而言，這種方式效果更好，在裂縫非?？拷鼏卧吔鐣r(shí)不會(huì)產(chǎn)生后者容易出現(xiàn)的數(shù)值誤差。但無論是常規(guī)有限元法在裂尖布置細(xì)密網(wǎng)格的方式，還是擴(kuò)展有限元法引入裂尖漸近位移場(chǎng)的方式，其主要目的都是為了提高裂尖附近的求解精度，從而提高應(yīng)力強(qiáng)度因子的計(jì)算精度，這些方法都有改進(jìn)的余地。即使是目前認(rèn)為最適合于裂縫擴(kuò)展分析的擴(kuò)展有限元法，由于其只是使用了裂尖漸近位移場(chǎng)的部分特征函數(shù)，嚴(yán)格地說，還不能構(gòu)成對(duì)裂尖位移場(chǎng)的最佳逼近。文獻(xiàn)[12]提出的新型數(shù)值流形方法，實(shí)現(xiàn)了裂紋尖端的解析解與其周邊數(shù)值解聯(lián)合應(yīng)用以求解應(yīng)力強(qiáng)度因子，能夠采用裂尖真實(shí)位移場(chǎng)的最佳逼近，并直接得到應(yīng)力強(qiáng)度因子，計(jì)算精度高，但僅限于平面問題的Ⅰ型和Ⅱ型裂紋。

在上述研究的基礎(chǔ)上，沿用解析解與數(shù)值解聯(lián)合應(yīng)用的思路，在裂紋尖端直接引入Williams解析解作為數(shù)學(xué)覆蓋，應(yīng)用裂紋尖端的解析解與周邊數(shù)值解的三維流形覆蓋聯(lián)系技術(shù)，可直接得出裂紋尖端的三維應(yīng)力強(qiáng)度因子，精度高且計(jì)算收斂快。

1 裂紋尖端位移場(chǎng)的Williams級(jí)數(shù)解

如圖1所示的裂紋尖端區(qū)域，Williams位移級(jí)數(shù)解

圖1 裂紋尖端坐標(biāo)系

(1)

(2)

2 裂紋尖端位移場(chǎng)的解析解覆蓋

在六面體網(wǎng)格內(nèi)，將式(1)寫成矩陣形式。

對(duì)于i=0:

(3)

對(duì)于i≥1:

(4)

將二維數(shù)值覆蓋強(qiáng)制約束方法[13]推至三維，形成獨(dú)立的三維解析覆蓋，權(quán)函數(shù)wj均為結(jié)點(diǎn)1的權(quán)函數(shù)w1。其中，m為覆蓋函數(shù)項(xiàng)數(shù)，下同。

對(duì)于第i項(xiàng)(i≥1)，其應(yīng)變子矩陣為

(5)

(6)

將各項(xiàng)偏導(dǎo)數(shù)展開

(7)

式中：

另外，

剛度子矩陣

(8)

式中：D為彈性矩陣；V為流形元體積。

3 裂紋周邊網(wǎng)格的解析解與數(shù)值解覆蓋

在裂紋周邊的各網(wǎng)格內(nèi)，位移統(tǒng)一表示為式(9)，其中i≥1。

(9)

式中：dnkl、enkl、fnkl為待求的多項(xiàng)式系數(shù)；n、k、l為多項(xiàng)式階次(p為多項(xiàng)式的最高階次)；wj1、wj2為相應(yīng)于各結(jié)點(diǎn)的權(quán)函數(shù)；J和L的個(gè)數(shù)根據(jù)網(wǎng)格而定，但保持J+L=8。顯然，取L=0則退化到解析解覆蓋式(4)。

如圖2所示(圖中窄條為部分重疊區(qū)域，長(zhǎng)方體區(qū)域?yàn)楠?dú)立覆蓋，陰影為裂紋)，采用強(qiáng)制約束的方法，將網(wǎng)格結(jié)點(diǎn)28、40、39、33、34、45、46約束到結(jié)點(diǎn)27，即令

圖2 六面體數(shù)學(xué)網(wǎng)格中的解析解覆蓋和數(shù)值解覆蓋

對(duì)于周邊的其他網(wǎng)格，比如：網(wǎng)格1-2-14-13-7-8-20-19中的結(jié)點(diǎn)均為數(shù)值解的獨(dú)立覆蓋；解析覆蓋的相鄰網(wǎng)格39-40-52-51-45-46-58-57中的結(jié)點(diǎn)39、40、45、46為解析結(jié)點(diǎn)，其他4個(gè)結(jié)點(diǎn)為數(shù)值結(jié)點(diǎn)；而裂紋自由端所在網(wǎng)格25-26-38-37-31-32-44-43在豎直方向采用了兩個(gè)獨(dú)立覆蓋，實(shí)現(xiàn)裂紋兩邊的獨(dú)立運(yùn)動(dòng)。

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

(10)

(11)

令f(x,y,z)=xn-k-lylzk，則式(11)可寫為

(12)

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

(13)

由式(10)和式(12)可見，剛度矩陣積分中具有非多項(xiàng)式的函數(shù)，基于多形式函數(shù)的單純形積分公式無法應(yīng)用，必須采用數(shù)值積分。因此，采用四面體區(qū)域的Hammer積分[14]，將流形元邊界上的三角片與其形心相連形成四面體，然后再進(jìn)行積分計(jì)算。

4 算例分析

以含邊界裂紋的無限長(zhǎng)柱體為例，考慮兩種不同類型荷載。

4.1 兩端受均布拉力

圖3 兩端受均布拉力的無限長(zhǎng)柱體內(nèi)的邊界裂縫

整體的流形元網(wǎng)格如圖4所示，共劃分6個(gè)獨(dú)立覆蓋和9個(gè)部分重疊覆蓋(窄條區(qū)域)，每個(gè)獨(dú)立覆蓋的大小基本相同。裂紋所在的獨(dú)立覆蓋區(qū)域大小為0.65 m×0.2 m×1 m(長(zhǎng)×寬×厚)，柱體底面施加法向約束。

圖4 部分重疊覆蓋的流形元網(wǎng)格

應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果見表1，考慮了裂紋尖端所在的獨(dú)立覆蓋取Williams級(jí)數(shù)的不同項(xiàng)數(shù)，獨(dú)立覆蓋周邊的數(shù)值解覆蓋可取多項(xiàng)式的不同階數(shù)。

表1 應(yīng)力強(qiáng)度因子

隨著覆蓋函數(shù)階數(shù)的增加，三維應(yīng)力強(qiáng)度因子計(jì)算值基本接近理論值。Williams級(jí)數(shù)的階數(shù)以及周邊數(shù)值覆蓋階數(shù)對(duì)計(jì)算值有較大影響。

1)Williams級(jí)數(shù)的階數(shù)(m)影響最大。當(dāng)m≤3，則KI與理論值相差較大；當(dāng)m≥4時(shí)，KI接近理論值。

2)周邊數(shù)值覆蓋階數(shù)升高有利于提高解的精度。當(dāng)周邊數(shù)值解取1階時(shí)，KI的計(jì)算值普遍小于取2階的情況，當(dāng)m≥4時(shí)，KI與理論值接近，但仍有差距，僅當(dāng)m取為7時(shí)才達(dá)到1.42，與理論值1.45最為接近。而當(dāng)周邊數(shù)值覆蓋取2階，m≥7時(shí)，計(jì)算值與理論值基本一致。

前期平面問題研究表明，在大的覆蓋中單純依靠提高覆蓋函數(shù)階次的方法往往會(huì)帶來計(jì)算結(jié)果的振蕩跳躍。反之，如果采用較小的覆蓋而用相對(duì)簡(jiǎn)單的低階多項(xiàng)式，則可以更好地逼近實(shí)際復(fù)雜的分布情況。表2的計(jì)算結(jié)果也說明，三維問題中，基于大覆蓋，僅僅采用提高級(jí)數(shù)解及相鄰覆蓋函數(shù)的階數(shù)的做法來提高計(jì)算精度，計(jì)算也表現(xiàn)出一定的不穩(wěn)定性，要取得滿意的計(jì)算精度，裂紋尖端及其周邊的覆蓋函數(shù)的階數(shù)不小于7階。

考慮到裂紋所在的獨(dú)立覆蓋區(qū)域較大，因此，將裂紋所在的獨(dú)立覆蓋進(jìn)行局部加密，將裂紋尖端覆蓋分別加密1倍及2倍，采用局部覆蓋加密技術(shù)[16]，重新計(jì)算應(yīng)力強(qiáng)度因子，如表2所示。

表2 應(yīng)力強(qiáng)度因子(n=2)

表2結(jié)果表明，當(dāng)裂紋尖端獨(dú)立覆蓋加密1倍后、Williams級(jí)數(shù)的階數(shù)≥4或當(dāng)裂紋尖端獨(dú)立覆蓋加密2倍后，Williams級(jí)數(shù)的階數(shù)≥3時(shí)，KI與理論解十分接近，且隨著Williams級(jí)數(shù)階數(shù)的提高，計(jì)算結(jié)果趨于穩(wěn)定。

4.2 兩端受剪切荷載作用

圖5 兩端受剪力的無限長(zhǎng)柱體內(nèi)的邊界裂縫

流形元網(wǎng)格如圖6所示，共劃分18個(gè)獨(dú)立覆蓋(方塊區(qū)域)和25個(gè)部分重疊覆蓋(窄條區(qū)域)，每個(gè)獨(dú)立覆蓋的大小基本相同。裂紋所在的獨(dú)立覆蓋區(qū)域大小為0.2 m×0.2 m×1 m(長(zhǎng)×寬×厚)，柱體底面施加法向約束。應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果見表3。

圖6 部分重疊覆蓋的流形元網(wǎng)格

表3 應(yīng)力強(qiáng)度因子

表3計(jì)算結(jié)果表明：采用圖6所示計(jì)算網(wǎng)格，當(dāng)m≥7時(shí)，周邊數(shù)值覆蓋階數(shù)取2、3的多項(xiàng)式階數(shù)時(shí)，計(jì)算結(jié)果與理論值比較符合。當(dāng)m≤7，計(jì)算結(jié)果與理論值有一定差別，局部數(shù)值出現(xiàn)跳躍，表明裂紋附近的網(wǎng)格還沒有達(dá)到足夠的密度。當(dāng)網(wǎng)格加密一倍后，周邊數(shù)值覆蓋階數(shù)均取2階，當(dāng)m≥7時(shí)，計(jì)算值與理論值十分接近，且隨著階數(shù)的提高，計(jì)算結(jié)果趨于穩(wěn)定。

以上算例驗(yàn)證了三維裂縫計(jì)算公式和程序的正確性，表明裂紋尖端解析解覆蓋和周邊數(shù)值解覆蓋聯(lián)合應(yīng)用求解三維線彈性斷裂力學(xué)問題可行。與常規(guī)有限元方法相比，無需在裂紋尖端布置細(xì)密的網(wǎng)格，計(jì)算精度高，收斂相對(duì)較快。

裂紋尖端獨(dú)立覆蓋的密度、解析覆蓋的級(jí)數(shù)以及相鄰數(shù)值覆蓋的階數(shù)是影響應(yīng)力強(qiáng)度因子計(jì)算精度的重要因素，但在保證獨(dú)立覆蓋有一定密度的情況下，提高與獨(dú)立覆蓋相鄰數(shù)值覆蓋的階數(shù)可以得到應(yīng)力強(qiáng)度因子的精確解。

裂紋尖端獨(dú)立覆蓋的合理布置對(duì)應(yīng)力強(qiáng)度因子的計(jì)算精度及穩(wěn)定性有一定的影響，因此，下一步要重點(diǎn)研究裂紋尖端附近的覆蓋自動(dòng)布置及密度問題，以保證方法的收斂性，便于開展三維裂縫擴(kuò)展的動(dòng)態(tài)模擬研究。

5 結(jié)論

將裂紋尖端解析解覆蓋和周邊數(shù)值解覆蓋聯(lián)合應(yīng)用，分析三維線彈性斷裂力學(xué)問題，得到以下主要結(jié)論：

1)在包含裂紋尖端的解析覆蓋中，應(yīng)用裂紋尖端附近的Williams位移解析解作為覆蓋函數(shù)，并采用高階多項(xiàng)式三維覆蓋函數(shù)與解析覆蓋的條形連接技術(shù)，實(shí)現(xiàn)了在解析覆蓋中直接求得裂紋尖端的三維應(yīng)力強(qiáng)度因子。

2)典型的張開型和撕開型的裂紋算例表明，應(yīng)力強(qiáng)度因子的計(jì)算精度較高。鑒于三維裂縫擴(kuò)展問題的復(fù)雜性，裂紋尖端周邊數(shù)值覆蓋階數(shù)以及獨(dú)立覆蓋網(wǎng)格密度對(duì)應(yīng)力強(qiáng)度因子計(jì)算精度的影響較二維問題更大。因此，協(xié)調(diào)獨(dú)立覆蓋密度、階數(shù)與周邊三維數(shù)值覆蓋階數(shù)的關(guān)系，來保證高精度求解收斂性的快速、穩(wěn)定是下一步研究的重點(diǎn)。

考慮到解析級(jí)數(shù)是裂尖附近真實(shí)位移場(chǎng)的最佳逼近，相比其他方法而言，可以認(rèn)為該方法在應(yīng)力強(qiáng)度因子求解方面逼近效果更好、收斂更快，同時(shí)，由于網(wǎng)格布置根據(jù)不同區(qū)域的精度要求，只在裂尖附近進(jìn)行覆蓋加密，因此，相比采用均勻網(wǎng)格的擴(kuò)展有限元而言，計(jì)算效率將有所提高，可以實(shí)現(xiàn)大規(guī)模計(jì)算。另外，應(yīng)力強(qiáng)度因子SIF本身就是裂尖解析級(jí)數(shù)的未知數(shù)，在求解系統(tǒng)方程組時(shí)一并得到，而不需要像其它方法那樣通過所謂的“直接”法或“間接”法來推求，不僅方便，而且不會(huì)引入額外誤差，這也是該方法的優(yōu)勢(shì)所在。

該方法可以同時(shí)求解Ⅰ型、Ⅱ型、Ⅲ型(撕開型)裂紋的應(yīng)力強(qiáng)度因子，應(yīng)用復(fù)合型裂紋擴(kuò)展準(zhǔn)則就可以判斷其是否繼續(xù)開裂，因此，該方法在三維裂縫擴(kuò)展的動(dòng)態(tài)模擬方面極具應(yīng)用前景。