基于3D-VCCT的多裂紋輸氣管道斷裂分析方法

姚安林 赫文博 徐濤龍 蔣宏業 谷東方

1.西南石油大學石油與天然氣工程學院 2.油氣消防四川省重點實驗室

0 引言

油氣管道在制造、安裝或服役過程中不可避免地會產生裂紋缺陷，這些缺陷有時是以多裂紋的形式存在的[1]。GB/T 19624—2004《在用含缺陷壓力容器安全評定》在一定條件下采用包絡合并準則將相鄰裂紋合并成一個裂紋再進行安全評定，而相關研究[2]認為包絡合并準則不能準確反映裂紋間的干涉影響，從而使得評定結果偏于保守。

目前，國內外學者對管道裂紋已開展了一定研究。Kumar等[3]以管道周向裂紋為研究對象，開展管道裂紋尖端應力強度因子（K）的研究；Jacquemin等[4]針對管道焊縫錯位裂紋，分析了裂紋深度和焊縫錯位值對斷裂參量J積分的影響；Arafah等[5]以厚壁管道為研究對象，開展雙軸拉伸彎曲載荷試驗，證明了預測臨界載荷和試驗所得載荷的吻合性。宋漢成等[6]借鑒管道軸向裂紋和環向裂紋的評估方法，提出了基于應力投影的螺旋裂紋評估方法；白永強等[7]通過對輸氣管道裂紋擴展因素進行分析，發現管道內壓對裂紋驅動力的影響最大，但以上研究均是基于單裂紋缺陷。而相對于單裂紋缺陷，管道多裂紋缺陷的建模難度更大、干涉機理更復雜，目前僅有少數學者進行了研究，何雪等[8]將含雙軸向裂紋缺陷的管段簡化為二維模型，討論了裂紋夾角對應力強度因子的影響，但該研究工作未對簡化模型的合理性作出解釋。綜上所述，針對輸氣管道多裂紋干涉效應的研究相對較少且不成熟，因此有必要繼續開展此類研究，以期提高含缺陷輸氣管道安全評估的準確性。

以斷裂力學為基礎的失效評估圖（Failure Assessment Diagram，FAD）為含裂紋缺陷構件的安全評定提供了新的途徑，其中K是運用FAD評估裂紋缺陷的重要斷裂參量之一[9]。鑒于此，首先引入奇異單元法和3D-VCCT的理論模型，給出K的計算方法；其次，在對比奇異單元法和3D-VCCT的基礎上，優選3D-VCCT建立輸氣管道多裂紋干涉模型；最后，通過分析裂紋相互干涉作用因子（η）的變化，研究附屬裂紋對主裂紋的干涉效應，從而為后續建立更精準的多裂紋安全評定模型提供技術指導和參考數據。

1 基本理論

對于承受內壓的輸氣管道而言，由環向應力引起的張開型（Ⅰ型）斷裂相對較為嚴重，故以軸向半橢圓裂紋為研究對象。目前K的計算方法主要有解析法、試驗法和數值法，其中解析法只能求解特定幾何形狀和荷載的單裂紋問題，試驗法雖可以求解K，但成本較高[10]。相比于解析法和試驗法，數值法在斷裂問題的求解過程中受到越來越多學者的青睞，其中3D-VCCT和奇異單元法是兩種重要的方法。

1.1 3D-VCCT理論模型

針對一個厚度為B的平板線狀裂紋，Irwin[11]提出能量釋放率（G）的概念，其定義為裂紋擴展（Δc）產生新裂紋面所需要的能量。對于Ⅰ型裂紋，其表達式為：

式中GI表示Ⅰ型裂紋的能量釋放率，N/mm；W表示外力做功，N·mm；B表示厚度，mm；Δc表示裂紋擴展量，mm；σ表示應力，N；u表示豎向（Y軸）距離，mm；r表示距裂紋尖端的距離，mm。

式（1）涉及沿著閉合裂紋線對應力的積分，而在有限元分析中虛擬裂紋線上應力做功等于節點力在節點位移上做功，Ⅰ型能量釋放率可以通過結點力與節點位移計算得到，其表達式如下：

式中u1,2表示節點1與節點2在豎向的位移變化，mm；Fy5表示節點5的節點力，N。

相對于線狀裂紋，面狀裂紋[12]由裂紋前緣而不只是一個點（線狀裂紋的裂紋尖端）來描述，因此Ⅰ型能量釋放率的計算相對更加復雜。有限元分析中面狀裂紋八節點單元的能量釋放率計算公式如下[13]：

式中 表示i單元的Ⅰ型能量釋放率，N/mm；wi表示i單元的裂紋前緣長度，mm；F1、F2分別表示節點力，N；u1、u2分別表示節點的位移變化量，mm；C1、C2分別表示系數；wi-1表示i-1單元的裂紋前緣長度，mm；wi+1表示i＋1單元的裂紋前緣長度，mm。

對于平面應變問題，應力強度因子與能量釋放率存在如下關系：

式中E表示彈性模量，MPa；μ表示泊松比； 表示i單元的Ⅰ型應力強度因子，MPa·mm1/2。

1.2 奇異單元法理論模型

Henshell和Shaw[14]于20世紀70年代中期提出了奇異單元理論，奇異單元是繞裂紋尖端附近構造的一種特殊單元，可反映裂紋尖端應力應變的奇異性。對于平面應變狀態，Shih等[15]證明了裂紋尖端Ⅰ型應力強度因子可通過下式求得：

式中r1/4表示1/4節點到裂紋尖端的距離，mm；v1/4表示1/4節點處裂紋的張開位移，mm。奇異單元的節點位置如圖1所示。

圖1 奇異單元的節點位置圖

圖2 輸氣管道裂紋幾何模型圖

2 管體三維半橢圓裂紋模型

2.1 幾何模型及參數設定

輸氣管道裂紋幾何模型如圖2所示，選取管長（L）為2 000 mm，壁厚（t）為17.5 mm，內徑（R）為490.5 mm，均布內壓載荷（p）為10 MPa的X70管線鋼為研究對象[16]，其力學參數為：泊松比0.3；密度7 900 kg/m3；彈性模量210 GPa。圖2-b給出了平行共線、平行共軸和平行偏置3種位置關系的多裂紋幾何模型。設定主裂紋長半軸為c，短半軸為b，裂紋前緣與內壁夾角為θ，θ＝0°為近表面點，θ=90°為最深點，θ=180°為遠表面點，遠表面點與近表面點統稱為表面點，附屬裂紋長半軸為ci，短半軸為bi，主裂紋與附屬裂紋最深點的水平間距為Li，最深點的豎向間距為βi，i=1表示平行共線位置，i=2表示平行共軸位置，i=3表示平行偏置位置（圖2-b）。

2.2 有限元模型建立

基于ABAQUS有限元軟件，運用奇異單元法和3D-VCCT建立管道三維有限元模型，模型兩端施加固定約束，為消除邊界條件對裂紋前緣應力分布的影響，半橢圓裂紋預制在管道的中間位置，并在管道內壁和半橢圓裂紋面分別施加內壓載荷10 MPa。值得注意的是，使用有限元軟件計算K需要對裂紋面做出合適的幾何切分以便劃分出特定的裂紋網格，考慮到奇異單元法中與裂紋前緣相接觸的單元需采用奇異單元，故裂紋面切分、掃略（圖3-a），而3D-VCCT無需考慮裂紋單元的奇異性，裂紋面切分如圖3-b所示。此外，有研究[17]提到裂紋網格尺寸影響求解應力強度因子的精度，為此對裂紋前緣兩側的區域分別進行網格加密處理，網格劃分均采用八節點六面體全積分單元，對應的有限元單元類型為C3D8，奇異單元法和3D-VCCT法裂紋前緣（黃色線條）兩側各加密6層的網格模型如圖4所示。

3 數值模擬與結果分析

3.1 有限元模型驗證

根據《應力強度因子手冊》[18]，對于內壁含一條軸向半橢圓裂紋的管體，其裂紋前緣應力強度因子（KⅠ）可表示為：

式中p表示管道內壓，MPa；R表示管道內徑，mm；t表示壁厚，mm；θ表示裂紋前緣某點與內壁夾角，（°）；c表示裂紋長半軸，mm；b表示裂紋短半軸，mm。

圖3 三維輸氣管道有限元模型圖

圖4 管道裂紋前緣網格模型圖

假設管道內壁有一軸向半橢圓裂紋，長半軸（c）為6.25 mm，短半軸（b）為5 mm，采用圖4中的奇異單元法和3D-VCCT建立有限元模型，將所得模擬結果和《應力強度因子手冊》計算結果對比，分析不同網格密度條件下兩種數值模擬方法的精度差異。定義不同網格密度下數值模擬結果與理論結果的相對誤差，可由下式得到：

式中f有限元值表示3D-VCCT或奇異單元法數值模擬結果；f理論值表示《應力強度因子手冊》計算結果。

從模擬結果和理論計算結果發現，半橢圓裂紋前緣各點的應力強度因子關于裂紋最深點對稱，故僅分析裂紋前緣0°～90°各點應力強度因子的相對誤差。圖5為3D-VCCT法和奇異單元法在不同網格密度下的相對誤差，A表示3D-VCCT裂紋前緣兩側各加密層，B表示奇異單元法裂紋前緣兩側各加密層。

圖5 裂紋前緣不同角度下的相對誤差圖

對比奇異單元法與3D-VCCT法在不同角度下的相對誤差值可知，3D-VCCT的計算精度明顯高于奇異單元法，并且當裂紋前緣兩側各加密6層網格時3D-VCCT法的計算精度最佳，此時最大誤差為1.45%。圖6為裂紋前緣兩側各6層網格時3D-VCCT的Mises應力結果（圖6-a、b）和奇異單元法的Mises應力結果（圖6-c、d），從圖6可以看出，兩種方法都能真實反映裂紋前緣的應力集中現象。但進一步對比兩種有限元模型的計算效率（表1）發現，3D-VCCT的計算效率較高。因此綜合考慮裂紋面切分難易程度、有限元計算精度和計算效率，最終選擇高精度、高效率的3D-VCCT方法、裂紋前緣加密6層網格作為后續多裂紋干涉效應分析的基礎。

圖6 3D-VCCT法與奇異單元法的Mises應力對比圖

表1 3D-VCCT與奇異單元法的計算效率表

3.2 多裂紋干涉效應表征

為了定量表征附屬裂紋對主裂紋的干涉效應，這里引入裂紋相互干涉作用因子（η），其表達式為：

式中Km(θ, L, β)表示裂紋水平間距（L）和豎向間距（β）下主裂紋前緣在角度（θ）處的應力強度因子；Ks(θ)表示管道相同位置存在同尺寸單裂紋時裂紋前緣在角度（θ）處的應力強度因子。當η＞1時，表示附屬裂紋的存在增大了主裂紋前緣的應力強度因子，具有增強效應；反之，當η＜1時，則表示附屬裂紋對主裂紋具有減弱效應；當η=1時，表示附屬裂紋的存在對主裂紋不產生影響；當∣η-1∣＜1%時，可認為附屬裂紋對主裂紋的干涉效應忽略不計。

3.3 結果分析與討論

3.3.1 裂紋形狀對應力強度因子最大值點的影響

圖7給出了長半軸（c）為6.25 mm，短半軸（b）為5 mm半橢圓裂紋前緣（0°～90°）各點的應力強度因子。從圖7可以看出，應力強度因子在裂紋表面點和最深點存在極大值。為進一步探討裂紋形狀比對應力強度因子最大值點位置的影響，在圖3-b有限元模型的基礎上，保持裂紋短半軸（b）為5 mm不變，計算=0.50、0.77、0.78、0.80、1.00和 1.50情況下的值（ 表示裂紋表面點的應力強度因子，表示裂紋最深點的應力強度因子），其結果如圖8所示。

圖7 裂紋前緣各點的應力強度因子圖

圖8 裂紋形狀比對應力強度因子最大值點位置的影響圖

3.3.2 裂紋尺寸對主裂紋的干涉效應

在獲取裂紋形狀比對應力強度因子最大值點位置的影響規律后，設定 ＜0.8并在此條件下建立3D-VCCT多裂紋干涉有限元模型，其中材料參數、載荷情況和邊界條件均保持不變，輸氣管道多裂紋網格干涉模型如圖9所示。

圖9 輸氣管道多裂紋網格模型圖

圖10給出了當裂紋處于平行共線位置關系時附屬裂紋尺寸變化對η的影響。由圖10-a可看出，η隨附屬裂紋長度的增加而增加，而隨主裂紋前緣各點角度的增加呈指數形式減小，當主裂紋前緣角度大于90°后，η值趨于穩定。由圖10-b可看出，η隨附屬裂紋深度的增加而增加，而隨主裂紋前緣角度的增大呈現先減少后略微增大的變化趨勢，主裂紋最深點的增強效應最弱。綜合圖10-a、10-b發現，當附屬裂紋存在時，主裂紋前緣各點的η均大于1，即主裂紋前緣各點受到了附屬裂紋的增強效應，且近表面點的增強效應最明顯，這是因為主裂紋近表面點與附屬裂紋相距最近，應力疊加最嚴重（圖11），故在對平行共線裂紋進行安全評定時需重點考慮近表面點。

圖12分別給出了當裂紋處于平行共軸位置關系時，主裂紋前緣各點η隨附屬裂紋長度和深度的變化情況。由圖12-a可看出，η在主裂紋最深點和表面點表現出不同的變化規律，即主裂紋最深點的η隨附屬裂紋長度的增加而增大，而主裂紋表面點的η卻隨附屬裂紋長度的增加而減少。

圖10 平行共線位置關系下附屬裂紋尺寸變化對相互干涉作用因子的影響圖

圖11 平行共線位置關系下主裂紋與附屬裂紋的應力干涉云圖與單裂紋應力對比云圖

圖12 平行共軸位置關系下附屬裂紋尺寸變化對相互干涉作用因子的影響圖

由圖12-b可知，當裂紋間的豎向間距較大時，附屬裂紋對主裂紋前緣各點的∣η-1∣＜1%，裂紋間的干涉影響可忽略不計；但隨著附屬裂紋深度的增加，裂紋豎向間距不斷減小，主裂紋表面點開始產生減弱效應，且減弱效應隨附屬裂紋深度的增加而增強；值得注意的是，當b2/t≥0.3時，主裂紋最深點開始產生明顯的增強效應。從圖12-a和圖12-b不難發現，對于平行共軸裂紋，裂紋之間的相互干涉產生增強效應或者減弱效應，當裂紋豎向間距小于一定數值后主裂紋表面點產生減弱效應，而最深點產生增強效應，并且干涉效應隨附屬裂紋深度或長度的增加而增強。圖13為平行共軸位置關系下主裂紋與附屬裂紋相互干涉的Mises應力云圖與單裂紋應力云圖的對比，從圖13可以看出，由于雙裂紋的相互干涉效應，主裂紋表面點的應力明顯降低，而最深點的Mises應力在增大，因此在進行安全評定時需重點考慮裂紋最深點。

3.3.3 裂紋間距變化對主裂紋的干涉影響

為揭示裂紋間距對相互干涉作用因子的影響，保持主裂紋和附屬裂紋尺寸不變，通過變化裂紋水平間距（L）來開展多裂紋相互干涉效應研究，其結果如圖14所示。

圖13 平行共軸位置關系下主裂紋與附屬裂紋的應力干涉云圖與單裂紋應力對比云圖

圖14 裂紋距離變化對主裂紋相互干涉作用因子的影響圖

從圖14-a可知，當裂紋處于平行共線位置關系時，主裂紋前緣各點的η隨著水平間距L1的增大而逐漸減少并趨近于1，說明隨著裂紋水平間距的增大，裂紋前緣的應力逐漸得到釋放，從而導致主裂紋受附屬裂紋的增強效應愈來愈弱。

從圖14-b可以看出，當裂紋處于平行偏置位置關系時，隨水平間距的增大，主裂紋前緣各點的η逐漸增大并趨近于1，這是因為隨裂紋間距的增大，主裂紋與附屬裂紋相互遠離，應力屏蔽效應減弱。綜合圖14-a和圖14-b發現，當裂紋水平間距L≥6c時，主裂紋前緣各點的∣η-1∣＜1%，此時在對多裂紋進行安全評定時可不考慮附屬裂紋對主裂紋的干涉影響。

4 結論

1）單一半橢圓裂紋應力強度因子（K）的最大值點與裂紋形狀比有關，隨著的增加，的值逐漸增大，而當＜0.8時，K的最大值點由裂紋表面點轉移至裂紋最深點。

2）對于平行共線裂紋而言，附屬裂紋的存在將使主裂紋前緣各點產生增強效應，且隨附屬裂紋深度或長度的變化，主裂紋近表面點的增強效應最為顯著。

3）對于平行共軸裂紋而言，附屬裂紋的存在不僅存在增強效應也存在減弱效應，當裂紋豎向間距較小時，主裂紋表面點表現為減弱效應而最深點表現為增強效應，且干涉效應隨附屬裂紋深度或長度的增加而增強。

4）當裂紋間的水平間距大于主裂紋長半軸的6倍時，平行共線裂紋與平行偏置裂紋的干涉效應微乎其微，此時在對多裂紋輸氣管道進行安全評定時可不考慮裂紋間的相互影響，直接將多裂紋簡化為單裂紋處理。