基于新型流形法的三維應力強度因子求解

祁勇峰，蘇海東，龔亞琦

(長江科學院 材料與結構研究所；水利部水工程安全和病害防治工程技術研究中心，武漢 430010)

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

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

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

1 裂紋尖端位移場的Williams級數(shù)解

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

圖1 裂紋尖端坐標系

(1)

(2)

2 裂紋尖端位移場的解析解覆蓋

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

對于i=0:

(3)

對于i≥1:

(4)

將二維數(shù)值覆蓋強制約束方法[13]推至三維，形成獨立的三維解析覆蓋，權函數(shù)wj均為結點1的權函數(shù)w1。其中，m為覆蓋函數(shù)項數(shù)，下同。

對于第i項(i≥1)，其應變子矩陣為

(5)

(6)

將各項偏導數(shù)展開

(7)

式中：

另外，

剛度子矩陣

(8)

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

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

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

(9)

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

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

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

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

對于解析覆蓋結點j處的第i項，其應變矩陣為

(10)

(11)

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

(12)

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

(13)

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

4 算例分析

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

4.1 兩端受均布拉力

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

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

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

應力強度因子的計算結果見表1，考慮了裂紋尖端所在的獨立覆蓋取Williams級數(shù)的不同項數(shù)，獨立覆蓋周邊的數(shù)值解覆蓋可取多項式的不同階數(shù)。

表1 應力強度因子

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

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

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

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

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

表2 應力強度因子(n=2)

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

4.2 兩端受剪切荷載作用

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

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

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

表3 應力強度因子

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

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

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

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

5 結論

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

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

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

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

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