納入裂尖位置實(shí)時(shí)預(yù)測的天然氣管道裂紋擴(kuò)展模擬方法

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

(1.中國石油大學(xué)(華東)儲(chǔ)運(yùn)與建筑工程學(xué)院，山東青島 266580；2.中國石油大學(xué)(華東)山東省油氣儲(chǔ)運(yùn)安全重點(diǎn)實(shí)驗(yàn)室，山東青島 266580)

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

1 迭代加載法

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

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

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

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

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

2.1 管道有限元模型

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

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

2.2 基于納入裂尖位置實(shí)時(shí)預(yù)測的迭代加載法載荷施加

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

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

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

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

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

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

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

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

圖5 試驗(yàn)與模擬結(jié)果對(duì)比Fig.5 Comparison of test and simulation results

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

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

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

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

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

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

4.1 迭代次數(shù)影響

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

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

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

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

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

4.2 氣體減壓模型

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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