甄 瑩，常 群，曹宇光，祖毅真，鈕瑞艷

(1.中國石油大學(華東)儲運與建筑工程學院，山東青島 266580；2.中國石油大學(華東)山東省油氣儲運安全重點實驗室，山東青島 266580)

延性斷裂作為制約高鋼級天然氣管道發(fā)展的瓶頸問題[1-3]，亟需對其開展深入研究。全尺寸爆破試驗作為最可靠的研究方法，因受其試驗費用、周期及難度等方面的限制而無法廣泛開展[4]，數(shù)值模擬方法成為協(xié)助其開展相關研究的有效手段。管道裂紋擴展過程較為復雜[5-11]，使用基于全尺寸爆破試驗數(shù)據(jù)擬合的氣體減壓模型結合成熟損傷斷裂模型的有限元技術可以簡化該過程中氣體減壓與結構變形間復雜的流固耦合[8,12-14]，在計算難度與效率方面具有較大的優(yōu)勢，但仍需解決幾個問題。一是選擇高效且準確的損傷斷裂模型來實現(xiàn)裂紋動態(tài)擴展過程的模擬,二是復雜加載的實現(xiàn)。隨著裂紋不斷擴展，以裂尖為分界的內壓分布變化必然引起邊界條件的連續(xù)變化，而重新定義每次增量的壓力分布是不切實際的[15- 16]，目前并未形成較好的處理方法。針對這些問題，筆者選用在斷裂韌性表征方面較為優(yōu)越的GTN模型，通過歸納大量全尺寸管道爆破試驗數(shù)據(jù)，提出管道裂紋擴展過程中裂尖位置隨加載時間分階段近線性變化的基本假設，根據(jù)該假設構建納入裂尖位置實時預測的迭代加載法；在此基礎上建立全尺寸管道有限元模型，結合氣體減壓模型及所提出的加載法對X80管道裂紋動態(tài)擴展過程進行模擬，并利用全尺寸爆破試驗數(shù)據(jù)對模型進行驗證。

1 迭代加載法

使用氣體減壓模型的管道裂紋擴展模擬需要將管內加載區(qū)域分為裂尖前和裂尖后兩個部分，隨著裂紋不斷擴展，以裂尖為分界的內壓分布變化必然引起邊界條件的連續(xù)變化。歸納了不同鋼級、不同尺寸管道全尺寸爆破試驗中裂尖位置隨時間的變化規(guī)律[16-19]如圖1所示。

圖1 全尺寸管道爆破試驗數(shù)據(jù)及其分段線性擬合Fig.1 Experimental data of full-scale pipeline burst test and its piecewise linear fitting

由圖1可知，不同管道裂尖位置均隨時間呈近線性變化，為使全尺寸管道爆破試驗數(shù)據(jù)集中的分布于線性擬合線附近，將整個裂紋擴展過程分為階段Ⅰ、階段Ⅱ及階段Ⅲ，對3個階段分別進行線性擬合(圖1)?？梢?，所有試驗數(shù)據(jù)基本位于擬合線上，擬合系數(shù)均大于0.99，說明在一定范圍內裂尖位置基本隨加載時間線性變化。因此提出管道裂紋動態(tài)擴展過程中裂尖位置隨時間分階段近似線性變化的基本假設。

根據(jù)上述假設，若管道裂紋擴展過程中擴展速度已知，則可近似得到裂尖位置隨時間的變化關系；然而，實際管道爆破試驗或管道斷裂事故的裂紋擴展速度難以預測。同時由于裂紋起裂及止裂階段線性規(guī)律并非嚴格滿足，若將管道全程裂紋擴展速度設為定值，則會產(chǎn)生較大誤差。因此提出一種納入裂尖位置實時預測的迭代加載法(定義為CTPI method)：將整個裂紋擴展過程離散為幾個階段，每個階段計算結束后，利用Python提取裂尖位置與時間數(shù)據(jù)并計算裂紋擴展速度，之后將裂尖位置表示為裂紋擴展速度與時間之函數(shù)關系，將此關系代入下一階段加載子程序中作為加載的邊界條件進行相應模擬。如此進行不同階段模擬的迭代，直到計算完成或裂紋充分擴展，則可近似實現(xiàn)每個增量步中邊界條件隨裂尖位置的實時變化。

2 天然氣管道裂紋動態(tài)擴展數(shù)值模擬

2.1 管道有限元模型

開展單軸拉伸與三點彎試驗，以獲得材料本構與GTN模型損傷參數(shù)，相關試驗與參數(shù)標定見文獻 [20]。之后根據(jù)試驗所用管道尺寸(外徑D為1 422 mm，壁厚t為27.7 mm)，考慮其對稱性，建立四分之一模型，引入與管道外徑尺寸相同的初始裂紋。為保證裂紋可以充分擴展，將管道軸向總長的一半設為6倍外徑?？紤]到GTN模型嚴重的網(wǎng)格依賴性，將損傷斷裂過程區(qū)管壁厚度方向劃分12層單元，單元高度設為0.25 mm，采用過渡網(wǎng)格技術對管道網(wǎng)格進行劃分以實現(xiàn)單元密度沿軸向及環(huán)向的逐漸遞減，以此縮短計算時間。

有限元模型軸向采用了兩種不同尺寸的單元(軸向單元長度Le=0.5 mm和Le=1 mm)來研究網(wǎng)格尺寸效應。圖2為Le=0.5 mm時對稱顯示為二分之一管道的有限元模型。將已標定的GTN損傷演化參數(shù)賦予管道斷裂過程區(qū)，其余位置設置傳統(tǒng)彈塑性材料參數(shù)。

2.2 基于納入裂尖位置實時預測的迭代加載法載荷施加

Nonn等[19]提供了與本文所研究X80管道材料性能相似、尺寸相同的全尺寸管道爆破試驗數(shù)據(jù)(圖1(d))，依據(jù)該試驗工況對管道進行加載，并以試驗數(shù)據(jù)驗證模型可靠性。考慮模型的對稱性，在管道四分之一模型靠近初始裂紋一端施加z方向的對稱約束，軸向施加y方向的對稱約束，同時管道另一端施加等效拉力模擬端蓋效應。

圖2 全尺寸管道有限元模型Fig.2 Finite element model of full-scale pipeline

圖3 二維指數(shù)型氣體減壓模型Fig.3 2D exponential gas decompression model

利用提出的納入裂尖位置實時預測的迭代加載法對管道內部進行加載，將管道變形路徑劃分為多個階段，在初始階段以原始裂尖位置作為加載的邊界，借助Fortran語言編寫VDLOAD子程序，將加載區(qū)域分為裂尖前與裂尖后兩個區(qū)域。裂尖前根據(jù)試驗管道設計系數(shù)0.72及模擬所用材料本構參數(shù)，確定初始壓力為16.38 MPa，加載過程中該數(shù)值在20 ms內線性下降至初始水平p0的40%；裂尖后選用忽略壓力沿環(huán)向變化的二維指數(shù)型氣體減壓模型對管道進行加載(圖3)。初始模擬完成后，利用Python提取裂紋擴展速度，并將裂尖位置表示為裂紋擴展速度與加載時間之函數(shù)關系，將此函數(shù)代入下一階段VDLOAD加載子程序中，以此預測新的裂尖位置并將其作為裂尖前后兩個加載區(qū)域的分界點。由于所代入的裂尖位置函數(shù)與時間相關，因此加載過程中邊界條件實時更新，如此進行不同階段模擬的迭代，直到計算完成或裂紋充分擴展。上述管道的邊界條件及載荷施加的示意圖如圖4所示。

圖4 管道模型邊界條件及載荷示意圖Fig.4 Schematic diagram of boundary conditions and load for pipeline model

3 天然氣管道裂紋動態(tài)擴展數(shù)值模擬結果

將管道裂尖位置和速度歷史的數(shù)值結果與X80管道全尺寸爆破試驗[19]的試驗結果進行對比，如圖5所示。可見，借助本文提出的迭代加載法，兩種單元有限元模型所得管道裂尖位置隨時間的變化基本吻合，與試驗結果趨勢一致，最大誤差約為10%，可以認為滿足工程需求。

由圖5可以看出，隨時間增加，兩模型裂紋擴展長度增加均變緩，裂紋擴展速度均隨之出現(xiàn)下降趨勢，即管道趨于止裂。Le為0.5和1 mm的模型結果差別較小，證明在當前網(wǎng)格劃分策略下模擬結果對網(wǎng)格不敏感，模型結果與試驗結果趨勢一致，說明數(shù)值模型中假定的氣體減壓與實際減壓行為沒有顯著差異，也進一步證明了本文中所提出的納入裂尖位置實時預測的迭代加載法可以實現(xiàn)邊界條件的實時移動，近似實現(xiàn)了氣體減壓、管道變形與裂紋擴展的多場耦合。采用Le為0.5 mm的模型進行后續(xù)分析。

圖5 試驗與模擬結果對比Fig.5 Comparison of test and simulation results

為了進一步說明納入裂尖位置實時預測的迭代加載法(CTPI method)的必要性，將該加載過程進行簡化：同樣采用分步迭代的加載方法，直接將上一步計算完成時的裂尖位置作為下一步加載的分界點，定義該方法為簡化的迭代加載法(S-CTPI method)；同時考慮氣體無減壓的情況，即給管道施加均布壓力，不考慮壓力隨時間與位置的變化，定義該方法為均布加載法(NGD method)。采用這3種加載方法對管道裂紋動態(tài)擴展過程進行模擬，將每一階段預測并作為加載邊界條件的裂尖位置與模擬所得實際裂尖位置隨時間的變化關系匯總于圖6。

為便于比較，將圖6曲線劃分為3個階段?？梢姡谝浑A段裂紋擴展速度逐漸加快，但時間較短，CTPI method與S-CTPI method數(shù)據(jù)吻合較好，但均布加載法NGD method所得裂紋擴展長度均大于上述兩種方法。第二階段裂紋位置變化與時間基本呈線性關系，該階段裂紋擴展速度基本穩(wěn)定，因此CTPI method與S-CTPI method均分兩次迭代。從圖6中可以看出，CTPI method中由前一階段裂尖擴展速度預測并代入后一階段壓力分布計算的裂尖位置與實際模擬所得裂尖位置吻合程度極高；盡管S-CTPI method未考慮裂尖位置的連續(xù)變化，但仍得到了連續(xù)的模擬結果，且結果與考慮了加載邊界實時變化的CTPI method基本吻合。不考慮氣體減壓的NGD method與兩種迭代加載法模型數(shù)據(jù)也較為吻合。第三階段較第二階段時間短，但由于裂紋擴展速度不再保持定值，該階段均分3次迭代。由圖6可知，CTPI method模擬所得該階段裂紋尖端位置變化曲線與穩(wěn)定擴展階段相比增速已明顯下降，對應于圖5中裂紋擴展速度下降階段，該階段預測并應用的裂尖位置與模擬所得裂尖位置仍然吻合，說明CTPI method可以有效模擬管道動態(tài)裂紋的止裂現(xiàn)象。S-CTPI method模擬所得裂尖位置仍保持與裂紋穩(wěn)定擴展階段一致的增速，裂紋擴展速度未下降，而未考慮氣體減壓的NGD method所得裂尖位置變化曲線斜率甚至呈現(xiàn)上升趨勢。

圖6 不同加載方式下裂尖位置隨時間變化Fig.6 Change of crack tip position with time under different loading methods

進一步計算3種加載方法所得裂紋擴展速度如圖7所示。由圖7可以看出，CTPI method所得裂紋擴展速度曲線最為平滑，而另外兩種方法曲線均存在突變。CTPI method所得裂紋擴展速度在裂紋穩(wěn)定擴展階段內基本保持定值，而后隨氣體減壓呈現(xiàn)下降趨勢；S-CTPI method模擬所得裂紋擴展速度均在達到最大值后又在小范圍內波動，未出現(xiàn)明顯下降趨勢；而未考慮氣體減壓的NGD method所得裂紋擴展速度呈現(xiàn)輕微上升趨勢，這是由于裂紋擴展長度足夠大，實際裂尖位置靠近管道末端，管道內壓和等效軸向拉力共同作用使得裂紋擴展速度增加。數(shù)值方面，CTPI method、S-CTPI method和NGD method所得裂紋擴展速度穩(wěn)態(tài)值或最大值依次增大。綜上所述，本文中所提出的CTPI method可以有效模擬管道動態(tài)裂紋的止裂及氣體減壓過程，而S-CTPI method和未考慮氣體減壓的NGD method無法準確描述裂紋擴展速度的下降，因此無法有效實現(xiàn)管道止裂過程的模擬。

圖7 不同加載方式下裂紋擴展速度隨時間變化Fig.7 Change of crack velocity with time under different loading methods

圖8 不同加載方式下管道塑性變形Fig.8 Plastic deformation of pipeline under different loading methods

模擬所得管道變形情況見圖8。由圖8可以看出，塑性變形主要集中于裂紋張開的襟翼上，在相同時刻，CTPI method、S-CTPI method及NGD method所得管道裂紋張開程度及裂紋附近塑性應變大于10%的區(qū)域依次增加。CTPI method所得變形圖與實際管道試驗更為接近，襟翼顯示了移動的裂紋尖端后連續(xù)的“波浪”狀變形；S-CTPI method所得管道裂紋張開程度較大，裂紋襟翼上殘余了大量塑性變形，這是由于模型加載子程序中使用的作為加載邊界的裂紋尖端位置未能實現(xiàn)實時移動，加載邊界位于實際裂尖后某位置處，使部分管道未能遵循氣體減壓模型的壓力遞減規(guī)律而長時間承受較大內壓；未考慮氣體減壓的NGD method模擬所得管道變形最大，裂尖后整個區(qū)域發(fā)生了完全屈服，管壁張開程度超過180°，當加載時間達到0.018 s時，管壁發(fā)生了嚴重畸變，已不具備參考價值。由此可知，若不考慮氣體減壓則無法有效實現(xiàn)管道裂紋擴展的有效模擬。

4 管道裂紋動態(tài)擴展模擬方法

4.1 迭代次數(shù)影響

不同迭代次數(shù)時模擬所得裂尖位置變化曲線見圖9。由圖9可以看出，在每個迭代步內，由前一階段預測并嵌入加載子程序作為加載分界點的裂尖位置與模擬所得實際裂尖位置之差由零逐漸增加，在該迭代步結束時達到最大。對于五次迭代情況，該差值整體較小，最大僅為0.15 m，證明加載與實際情況吻合較好；對于二次迭代情況，在高壓力水平下該差值最大為0.8 m，在低壓力水平下最大為0.23 m；而且裂尖位置預測值始終大于模擬值，導致裂尖處及裂尖后部分位置位于減壓區(qū)，使得所受壓力小于實際情況而開裂壓力減小，從而使裂紋擴展速度過早下降。

圖9 迭代次數(shù)對裂尖位置影響Fig.9 Effect of iteration number on crack tip position

圖10為不同迭代次數(shù)時由前一階段預測并嵌入加載子程序作為加載分界點的裂尖位置與模擬所得實際裂尖位置之差。由圖10可以看出，迭代次數(shù)越少，該差值最大值越大。圖11為迭代次數(shù)對裂紋擴展速度的影響。由圖11可以看出，迭代次數(shù)對初始階段和穩(wěn)定裂紋擴展階段的裂紋速度影響不大。對于第三階段，迭代次數(shù)越少，裂紋速度下降越明顯。對于本文模擬工況，當?shù)螖?shù)達到四次時，模擬結果趨于穩(wěn)定，且不隨迭代次數(shù)增加而變化。

圖10 迭代次數(shù)對裂尖位置預測值與模擬值之差影響Fig.10 Effect of iteration number on differences between predicted and simulated crack tip position

圖11 迭代次數(shù)對裂紋擴展速度影響Fig.11 Effect of iteration number on crack velocity

4.2 氣體減壓模型

氣體減壓模型為對現(xiàn)有為數(shù)不多的真實管道爆破試驗數(shù)據(jù)和前人數(shù)值計算結果的簡化與近似。使用的氣體減壓模型主要包括二維指數(shù)型(圖3)、二維線性(圖12)、三維指數(shù)型(圖13)3種，圖14為3種模型模擬所得裂尖位置與裂紋擴展速度。

圖12 二維線性氣體減壓模型Fig.12 2D linear gas decompression model

圖13 三維指數(shù)型氣體減壓模型Fig.13 3D exponential gas decompression model

由圖14可以看出，兩種指數(shù)型模型模擬所得結果在整個裂紋擴展階段表現(xiàn)出良好的一致性，曲線基本重合；二維線性模型所得數(shù)據(jù)在裂紋穩(wěn)定擴展階段中期之前與兩種指數(shù)型模型吻合較好，但隨著持續(xù)加載，其模擬所得裂尖位置仍保持近似線性增長的趨勢，裂紋擴展速度也并未下降。說明在管道動態(tài)裂紋擴展模擬過程中，3種模型不可以相互替代，在相同壓力水平及其他設置下，線性模型所得裂紋擴展速度要大于其他兩種模型，所得止裂長度也會更大，將其用于止裂設計研究可能會出現(xiàn)非保守的結果。

圖14 氣體減壓模型對裂尖位置和裂紋擴展速度影響Fig.14 Effect of gas decompression model on crack tip position and crack velocity

圖15為不同氣體減壓模型模擬所得管道塑性變形情況。由圖15可以看出，在當前壓力水平下，二維與三維指數(shù)型模型所得管道變形與塑性分布情況幾乎一致，二維線性模型得到了較大的裂紋開口，裂紋尖端周圍塑性應變大于10%的區(qū)域也明顯大于兩種指數(shù)型模型。3種氣體減壓模型均實現(xiàn)了裂紋尖端后“波浪”狀變形的模擬，其中二維線性模擬所得“波浪”變形更為密集，說明其軸向拉伸效應更顯著。二維與三維指數(shù)型模型所得裂紋擴展長度、擴展速度、管道變形及塑性應變分布情況均一致，因此二者相互替代，考慮模型形式更為簡單，推薦使用二維指數(shù)型氣體減壓模型對管道動態(tài)裂紋擴展進行模擬。

圖15 氣體減壓模型對管道塑性變形的影響Fig.15 Effect of gas decompression model on plastic deformation of pipeline

5 結 論

(1)納入裂尖位置實時預測的迭代加載法可有效實現(xiàn)裂紋擴展過程中加載邊界的實時移動，近似實現(xiàn)了氣體減壓、管道變形與裂紋擴展的多場耦合，可以有效反應管內氣體分布情況及其隨裂尖移動的實時更新;應用該方法可有效實現(xiàn)裂紋擴展及止裂模擬，所得管道斷裂參量、裂紋張開程度及裂紋尖端后典型的連續(xù)的“波浪”狀變形與實際管道試驗結果類似。

(2)納入裂尖位置實時預測的迭代加載法一定程度上依賴于迭代次數(shù)，在每個迭代步內，由前一步預測并嵌于加載子程序作為加載分界點的裂尖位置與模擬所得實際裂尖位置之差由零逐漸增加，在該迭代步結束時達到最大；理論上講，迭代次數(shù)越多，模擬越接近實際情況。

(3)二維與三維指數(shù)型氣體減壓模型模擬所得裂紋擴展長度、速度、管道變形均較為一致，二者可以相互替代，考慮二維模型形式較為簡潔，推薦將其用于后續(xù)天然氣管道裂紋動態(tài)擴展的模擬。