伍 穎 吳 鵬 張 媛 鄢小兵

1. 西南石油大學土木工程與測繪學院 2. 東方電氣集團東方鍋爐股份有限公司 3.中國石油西南油氣田公司川中油氣礦

0 引言

油氣管道在搬運和回填的過程中經常會發(fā)生管道碰撞或巖石等外物壓砸管道的情況，造成管道凹痕缺陷[1-8]。凹痕是指由于管壁永久塑性變形而使管道橫截面發(fā)生的總的變形，是造成長輸油氣管道失效的一個常見原因[9-10]。凹痕缺陷會改變管體局部應力、應變狀態(tài)[11]，并增加管道的失效敏感性，直接影響管道的局部應變、剩余強度和疲勞壽命[12]。因此，需要對油氣管道上的凹痕進行正確評估，以判斷是否對變形管段進行修復或更換[13-14]。

傳統(tǒng)的基于凹痕深度的評價標準，忽略了凹痕其他的幾何形狀特征，忽視了凹痕區(qū)域的局部應變和應力分布，會導致管道的安全狀況被高估或低估[15]。因此，Rosenfeld等[16]提出了基于應變的評價思想，美國機械工程師協(xié)會ASME B31.8-2003[17]將該思想納入其中，制訂了從應變的角度評估凹痕缺陷的準則，成為傳統(tǒng)深度準則的替代方案。為了獲得凹痕區(qū)域應變，ASME B31.8-2003[17]中給出了計算凹痕區(qū)域等效塑性應變的非強制表達式，隨后在ASME B31.8-2007[18]、ASME B31.8-2018[19]中進行了部分改進。Baker等[20]認為標準中沒有指明應變計算的假設條件和應變確定的依據，計算等效塑性應變的方法需要進一步改善。Lukasiewicz等[21]綜合采用數(shù)學方法和基于殼單元的有限元模型，針對管道特有的圓形薄殼結構進行了凹痕應變求解，但由于該模型考慮了管材的非線性本構關系，在涉及大變形問題時計算十分復雜。Tindall等[22]利用漏磁檢測數(shù)據來預測凹痕的軸向彎曲應變[23]，提出了將軸向彎曲應變作為衡量凹痕嚴重程度指標的預測模型，不過該模型只考慮了凹痕處的軸向彎曲應變，因此存在一定的局限性。馮慶善等[24]利用慣性測繪數(shù)據，提出了基于彎曲應變的管道凹痕建模計算方法，但該方法也只考慮了彎曲應變與凹痕深度的關系。Adrian等[25]利用管壁質點在發(fā)生變形時產生的位移進行應變求解，但作者假設管道發(fā)生凹痕變形時只存在徑向位移，這就使得計算結果缺乏可信度。Husain[26]使用概率統(tǒng)計方法推導了管道凹痕區(qū)域應力應變的數(shù)學解析式，但該方法公式過于復雜，不利于工程實際應用。筆者團隊也從凹痕輪廓插值和應變計算方面對凹痕評價進行了研究，通過彎曲應變使用關系表達式計算膜應變可以得到很好的結果，但是這只針對特定情況，不具有普遍性[27-28]。上述應變計算方法在一定程度上能夠得到誤差更小的凹痕區(qū)域的應變結果，但是存在假設條件過多、局限性高、計算復雜、不利于工程實際應用的問題。Chike等[29-32]針對無壓管道形成的約束凹痕，通過以位移函數(shù)來定義凹痕輪廓，提出了基于應變張量進行應變求解的方法，但針對有壓管道及非約束凹痕的適用性還有待進一步研究。綜上所述，目前國際上仍尚無精確的分析方法用于確定實際管道凹痕處的應變[33]，并且過去的研究都忽略了剪切應變的計算，導致凹痕處總應變的確定可能存在偏差。

筆者提出了一種考慮剪切應變的計算內壓作用下含約束型和非約束型凹痕管道應變的新方法。首先，采用有限元方法對管道上約束凹痕和非約束凹痕的產生過程進行了數(shù)值模擬，分析了管道內外壁等效塑性應變沿管道軸向、環(huán)向的分布規(guī)律，然后通過對凹痕變形區(qū)域進行三維插值生成定義凹痕表面的位移函數(shù)，最終采用以應變張量形式分析凹痕區(qū)域應變的解析模型，建立了考慮剪切應變的管道凹痕區(qū)域應變計算模型?？紤]剪切應變的計算模型能夠更為準確地反應管道凹痕區(qū)域的應變分布，能為評估凹痕缺陷的嚴重程度提供一定的參考。

1 有限元數(shù)值模型

1.1 管材參數(shù)

管材選用API 5L X80管線鋼，管材的特性參數(shù)如表1所示[34]。

表1 鋼管的材料性能表

1.2 有限元模型建立

采用ABAQUS軟件建立橢球形壓頭作用于管道的三維實體模型，橢球形壓頭如圖1所示。筆者重點在于提出一種考慮剪切應變的計算內壓作用下含約束型和非約束型凹痕管道應變的新解析方法，而不是研究凹痕尺寸對應變的影響規(guī)律，因此假設橢球形壓頭短軸a、c大小相等，長軸為b。壓頭模型采用離散剛體建立，管道采取八節(jié)點線性六面體實體單元C3D8R建立，采用各向同性硬化模型。以國際上常用的準靜態(tài)分析方法模擬凹痕的形成，準靜態(tài)模擬以分步或多步將載荷非常緩慢地加載到管道結構上，相對于動態(tài)模擬可以忽略慣性效應[35-36]。由于管道模型和荷載的對稱性，建立了1/4對稱模型。通過在網格分布中產生一些偏差來優(yōu)化數(shù)值模擬的計算時間，局部網格細化分配在凹痕區(qū)域，而粗網格應用于管道模型的其他區(qū)域。

圖1 橢球形壓頭圖

有限元模型及網格劃分情況如圖2所示。

圖2 有限元模型網格劃分圖

壓頭加載采用位移加載方式。對受壓工況下的含約束凹痕管道進行模擬時，通過合理簡化，主要載荷加載步驟為：①建立接觸；②施加位移荷載；③施加內壓。對于非約束凹痕，有限元分析中的加載步驟為：①建立接觸；②施加位移荷載；③卸載位移；④施加內壓。非約束凹痕由于壓頭的移除會導致凹痕區(qū)域的彈性變形恢復，再施加內壓會進一步導致凹痕區(qū)域的塑性變形得到部分恢復，因此凹痕深度會減小。其中，內壓采用均布壓力載荷的方式均勻施加在管道的內表面。約束凹痕及非約束凹痕的有限元加載步驟分別如圖3、4所示。

圖3 約束凹痕加載步驟圖

圖4 非約束凹痕加載步驟圖

1.3 有限元模型驗證

為了驗證所建有限元模型的可靠性，對Naghipour等[37-38]中的實驗進行了數(shù)值模擬分析，并將模擬結果與試驗值進行對比。試驗管材為X80鋼，試驗管道長度為290 mm，管徑為44 mm，壁厚為2 mm，壓頭為直徑20 mm的半球。筆者建立了與試驗管參數(shù)一致的有限元模型，試驗結果和數(shù)值模擬結果對比如表2所示。從表2可以看出，該數(shù)值模型與試驗結果具有較好的一致性。將試驗邊界條件誤差考慮在一定范圍內是合理的，因此可以認為所建立的數(shù)值模型是有效的。

表2 實驗結果和有限元結果比較表

1.4 有限元應變結果分析

1.4.1 非約束凹痕

建立橢球形壓頭按壓管道形成非約束凹痕的有限元模型，壓頭短軸（a=c）為120 mm，長軸（b）為400 mm，管道外徑（D）為1 016 mm，壁厚（t）為15.3 mm，管道長度為3 048 mm，設計壓力為10 MPa，管道幾何尺寸選自國內某在役長輸管道的實際數(shù)據。凹痕深度（d）范圍為1%D～12%D（即10.16～121.92 mm）。圖5、6分別顯示了不同凹痕深度下凹痕區(qū)域管道內外壁沿管道環(huán)向和軸向的等效塑性應變分布情況。

從圖5-a可知，沿管道環(huán)向，最大等效塑性應變點始終位于凹痕中心的管道內表面，但隨著距凹痕中心距離的增加，應變值的變化規(guī)律呈現(xiàn)較為頻繁的波動。從圖5-b可知，管道外壁的等效塑性應變隨著距凹痕中心距離的增加不斷減小，最大等效塑性應變位置同樣位于凹痕中心，但管道外壁的應變遠小于管道內壁的應變。

圖5 非約束凹痕沿管道環(huán)向的等效塑性應變變化圖

從圖6可知，沿管道軸向，不論管道內表面還是外表面，凹痕區(qū)域等效塑性應變值的變化趨勢相似，但外壁等效塑性應變小于管道內壁。對比圖5、6，由于非約束凹痕失去壓頭的限制，因此會在內壓作用下發(fā)生塑性回復，導致凹痕區(qū)域邊緣的應變發(fā)生波動。

圖6 非約束凹痕沿管道軸向的等效塑性應變變化圖

1.4.2 約束凹痕

建立橢球形壓頭按壓管道形成約束凹痕的有限元模型，有限元模型幾何尺寸不變。圖7、8分別顯示了不同凹痕深度下凹痕區(qū)域管道內外壁沿管道環(huán)向和軸向的等效塑性應變分布情況。

從圖7-a可知，等效塑性應變隨著距凹痕中心沿管道環(huán)向距離的增加而呈非線性變化，最大等效塑性應變點始終位于凹痕中心。從7-b可知，凹痕區(qū)域管道外壁的等效塑性應變隨著距凹痕中心的管道環(huán)向距離的增加而不斷減小直至0，最大等效塑性應變位置同樣位于凹痕中心，但管道外壁的應變遠小于管道內壁的應變。這可能是因為管道施加內壓后，由于受到壓頭的限制，管道外表面處于受壓狀態(tài)而管道內表面處于受拉狀態(tài)，因此導致管道內壁的塑性應變急劇增加。

圖7 約束凹痕沿管道環(huán)向的等效塑性應變變化圖

從圖8-a可知，凹痕區(qū)域管道內壁的等效塑性應變值隨軸向距離的增加而不斷變化，但與管道環(huán)向相比，管道軸向的塑性變形區(qū)域更小，在很短的一段距離內等效塑性應變急劇減小。從圖8-b可知，凹痕區(qū)域管道外壁的等效塑性應變值的變化規(guī)律與內壁相似，但由于外表面直接受到壓頭的限制作用，在內壓作用下更容易受到壓頭的影響而產生塑性變形。此外，對比同等條件下的非約束凹痕，含約束凹痕的管道產生了更大的應變。

圖8 約束凹痕沿管道軸向的等效塑性應變變化圖

2 解析模型

解析模型參照了Chike等[29-32]提出的基于應變張量進行應變求解的方法，在此基礎上建立了考慮剪切應變的管道凹痕區(qū)域應變計算模型。

2.1 凹痕輪廓插值

凹痕區(qū)域應變可以通過對在線檢測工具的數(shù)據點插值形成變形輪廓來估計。選擇3次B-樣條來插值形成凹痕區(qū)域連續(xù)且可微分的三維拓撲曲面。樣條曲線遵循由凹痕坐標定義的控制多邊形的形狀，并且不圍繞任何直線振蕩。相鄰樣條函數(shù)的一階導數(shù)在它們相遇的節(jié)點處相等，二階導數(shù)也設為零。在這些條件下，由該方程可以得到插值常數(shù)，以分段的方式應用在管道的凹痕變形區(qū)域，因此可獲得凹痕區(qū)域連續(xù)且可微分的輪廓表面。

以橢球形壓頭形成的6%D（60.96 mm）約束凹痕為例，通過直接從數(shù)值模型的節(jié)點提取矢量位置來獲得插值點，并將有限元節(jié)點坐標轉換到圓柱坐標系（R，θ，Z）中。用于內插凹痕表面的數(shù)據點對應于在圓周方向上具有64個傳感器的在線檢測工具的分辨率，并且沿著管道的縱軸以每15 mm的間隔獲得數(shù)據。從凹痕頂點所在的橫截面作為鏡像中心，分析了沿管道長度300 mm的變形。使用計算工具MATHEMATICA完成對數(shù)據點的插值，示例凹痕內表面的三維插值曲面如圖9所示。

圖9 6%D約束凹痕插值內表面圖

2.2 幾何變形分析

2.2.1 基本假設

用于插值凹痕輪廓的樣條函數(shù)定義了管道內表面的幾何變形，因此可以通過分析管道中表面的變形獲得管道在軸向、環(huán)向和徑向的位移。所做的基本假設為：①假設凹痕頂點處經歷的變形是純徑向的；②假設變形過程中管壁厚度沒有減少；③假設管道在環(huán)向上不可延伸。

2.2.2 位移分解

定義凹痕頂點所在的角位置在0 rad處，整個圓周方向的角位置范圍為[－π, π]。圖10顯示了變形管的中面橫截面的示意圖，變形管道中間表面半徑為Rm，R0代表截面未變形的初始中面半徑，即為管徑（D）與壁厚（t）之差的一半。若一個點從未變形輪廓中的m移動到變形輪廓中的m'，管道變形后質點的角度變化為φ，則質點所移動的位移與角變形（φ）和中表面的環(huán)向位移有關。

圖10 變形和未變形狀態(tài)下管道中表面的剖面圖

通過在線檢測得到的數(shù)據是變形管道內壁半徑（Ri）和每個檢測點的相應角度位置（θ）。管壁中間表面的半徑（Rm）可以被評估為管道內壁半徑和管道厚度方向分量的總和。使用下式評估該半徑：

管道中表面徑向位移（urm）為：

管道中表面的環(huán)向位移（uθm）為：

管道的整體環(huán)向變形與管壁的橢圓化和由管壁在環(huán)向軸線上的傾斜引起的管壁局部變形有關。作為變形的結果，變形后的矢量位置與原始位置之間的環(huán)向角度形成三角關系，圖11顯示了管道在環(huán)向（R，θ）平面的變形，其中er和eθ表示為管道在徑向和環(huán)向方向上的單位矢量，θθ表示變形管道的中間表面在環(huán)向方向上形成的偏移角度。

圖11 管道在環(huán)向（R, θ）平面的變形圖

管道內表面的位移是管壁厚度的函數(shù)，管內壁環(huán)向位移（uθ）為：

管內壁徑向位移（ur）為：

其中θθ的計算方法如下：

變形管道的中間表面在軸向方向上形成偏移角度（θz），在管道的中間表面處，軸向位移（uzm）為0。內表面軸向位移（uz）由下式計算：

其中管壁軸向變形角度（θz）為：

2.3 應變測量

變形輪廓坐標以徑向位置（R）、角位置（θ）和軸向位置（Z）表示。管道變形在柱坐標系中管道的全局位移場（u）的定義為：

式中ez表示管道軸向方向上的單位矢量。

式中r表示管道的可變半徑，mm，它是管道的假設中間表面半徑（R0）與管壁厚度之和。

根據拉格朗日應變張量的定義，可得到格林—拉格朗日應變張量（E）為：

上述應變公式包含非線性項，因而可以解釋與變形相關的大變形和旋轉。當進行線性應變或小應變測量時，假設變形體材料顆粒的位移無窮小，因此在變形過程中，物體的幾何形狀不會發(fā)生變化。小應變張量（e）的定義表示為：

假設所研究的應變方向均在主軸上，則應變矩陣中的對角線應變分量εrr、εθθ和εzz分別代表徑向、環(huán)向和軸向應變分量；非對角線應變分量εrθ、εθz和εrz分別代表剪切應變分量。

在軸向方向，當管道發(fā)生凹痕變形時，以管壁膜的延伸為特征在軸向上發(fā)生較大變形，因此采用非線性拉格朗日應變測量，能夠獲得由于管道的膜延伸引起的應變。軸向應變的數(shù)學表示為：

在環(huán)向上，由于假設管道在環(huán)向上不可延伸，即管道變形后環(huán)向的總長度仍然保持不變，因而在該軸上管道發(fā)生了小的旋轉和小變形，應使用線性小應變公式來評估環(huán)向應變，其數(shù)學表達式為：

由于凹痕區(qū)域發(fā)生較大的塑性變形，可以忽略彈性的體積變化，認為材料在塑性狀態(tài)時的體積不可壓縮，這意味著有效應變矩陣的跡的和是0，則徑向應變公式為：

將3個方向的應變分量和剪切應變進行分析組合，能夠得到管道凹痕段的等效總應變?yōu)椋?/p>

若不考慮剪切應變，用徑向、環(huán)向、軸向的應變作為主應變表示時，

3 結果分析

3.1 非約束凹痕

3.1.1 變形對比分析

以初始加載深度為6%D（60.96 mm）的非約束凹痕為例，通過解析模型求解整個凹痕管道區(qū)域的方向位移表達式，這些表達式的解由分析模型得到的等值線表示（圖12～14），得到的管內壁徑向位移、環(huán)向位移、軸向位移分布如圖12-a、13-a、14-a所示。數(shù)值模擬得到的凹痕管道的徑向位移、環(huán)向位移、軸向位移分布如圖12-b、13-b、14-b所示。

從圖12-a可以看出，解析模型預測的變形區(qū)域中心的徑向位移最大，最大徑向位移為25.00 mm，負號代表發(fā)生向內的變形。徑向位移由初始加載深度6%D（60.96 mm）減小到25.00 mm，這是因為壓頭移除后，凹痕區(qū)域發(fā)生回彈現(xiàn)象，導致凹痕深度減小。隨著距凹痕中心沿管道軸向、環(huán)向距離的增加，徑向位移不斷減小直至為0，與圖12-b中有限元預測的頂點最大徑向位移為25.85 mm的結果一致。

圖12 解析模型和有限元模型預測的徑向位移圖

從圖13-a可以看出，沿管道環(huán)向方向，凹痕區(qū)域肩部出現(xiàn)了明顯的橢圓化現(xiàn)象。沿管道環(huán)向方向，凹痕頂點另一側會出現(xiàn)對稱的環(huán)向位移分布，負號表明為管壁發(fā)生反向的延伸。從圖13-b可以看出，有限元計算得到的最大環(huán)向位移出現(xiàn)在凹痕區(qū)域的側部，且最大環(huán)向位移量為25.65 mm，遠大于解析模型。這可能是由于解析模型中采用了管道環(huán)向不可延長的假設，才導致管道的環(huán)向位移被低估。但是在凹痕評估中，影響凹痕缺陷嚴重程度的主要因素為徑向位移，環(huán)向位移對凹痕缺陷評估影響不大，所以本模型可以很好地對凹痕缺陷進行評估。

圖13 解析模型和有限元模型預測的環(huán)向位移圖

從圖14-a可以看出，解析模型得到的管道最大軸向位移出現(xiàn)在凹痕區(qū)域的側部，最大軸向位移為1.00 mm。變形區(qū)域的軸向位移沿凹痕頂點對稱分布，正負號代表管壁沿軸向延伸的方向。從圖14-b可以看出，數(shù)值模擬預測的最大軸向位移同樣位于凹痕區(qū)域的管道內壁，且最大軸向位移量為1.50 mm，有限元模型與解析模型的預測基本一致。

圖14 解析模型和有限元模型預測的軸向位移圖

3.1.2 應變對比分析

通過解析模型求解管道凹痕區(qū)域的應變張量，各分量的值在分析模型得到的等值線圖上表示（圖15～17）。解析模型預測得到的管道內表面的徑向應變分量、環(huán)向應變分量、軸向應變分量分別如圖15-a、16-a、17-a所示，有限元模擬得到的管道內壁徑向應變、環(huán)向應變、軸向應變分別如圖15-b、16-b、17-b所示。

圖15 解析模型和有限元模型中管內壁的徑向應變圖

從圖15-a可以看出，解析模型求得的管道內表面最大徑向應變位于凹痕中心，值為0.100，負值代表壓縮應變。隨著距凹痕中心沿管道軸向、環(huán)向距離的增加，徑向應變不斷減小直至為0。從圖15-b可以看出，有限元模擬得到的管道內表面最大徑向應變同樣位于凹痕中心，值為0.071，與解析模型預測結果存在較小差異。

從圖16-a可知，解析模型預測的管內壁最大環(huán)向應變的位置位于凹痕中心，值為0.010。圖16-b中觀察到有限元分析得到的管道內壁最大環(huán)向應變點偏于凹痕頂點，位于凹痕頂點所在管道橫截面上凹痕區(qū)域的側部，管道內壁最大環(huán)向應變?yōu)?.069，大于解析模型對內壁環(huán)向應變的預測。這可能是由于上述提到的解析模型中對圓管環(huán)向不可延長的假設使得管壁環(huán)向變形被低估，因此導致了解析模型預測了較小的環(huán)向變形而使得估計的環(huán)向應變值偏小。

圖16 解析模型和有限元模型中管內壁的環(huán)向應變圖

從圖17-a觀察到管內壁的凹痕頂點位置的軸向應變最大，值為0.090，沿管道環(huán)向和軸向隨著距離凹痕中心距離的增大，管道內表面各點的軸向應變呈不斷減小的趨勢。數(shù)值模擬得到的管道內壁最大軸向應變?yōu)?.067，同樣位于凹痕中心，但值略小于解析模型預測的軸向應變。

圖17 解析模型和有限元模型中管內壁的軸向應變圖

考慮不同等效塑性應變組合公式，由解析模型得到的等效塑性應變分布如圖18所示。圖18-a為不考慮剪切應變的等效塑性應變分布，圖18-b為考慮剪切應變的等效塑性應變分布。

圖18 等效塑性應變分布圖

由有限元模型得到的等效塑性應變分布如圖19所示。

圖19 有限元分析得到的管內壁等效塑性應變分布圖

從圖18-a可知，不考慮剪切應變時，解析模型預測的凹痕管道的等效塑性應變的最大值為0.100，位于凹痕中心。從圖18-b可知，考慮剪切應變時的凹痕管道等效塑性應變最大值也為0.100，同樣位于凹痕中心，但應變分布范圍有所不同。對比圖18-a、b，剪切應變對等效塑性應變的峰值影響較小，因此剪切應變在一般的等效塑性應變求解中可以忽略，但剪切應變會影響凹痕區(qū)域的應變分布。

從圖19可知，有限元分析得到的管道內壁最大等效塑性應變?yōu)?.108，位于凹痕區(qū)域中心，有限元結果與解析模型預測結果基本一致。因此可以認為，對于非約束凹痕的總體應變評估，解析模型能很好地預測凹痕管道的最大應變值和峰值位置。

3.2 約束凹痕

將上述步驟應用于本研究中的約束凹痕，得到了相似的變形和應變結果。由解析模型得到的深度為6%D（60.9 mm）的凹痕管道內表面等效塑性應變分布如圖20所示，由有限元模型得到的等效塑性應變分布如圖21所示。

圖20 等效塑性應變分布圖

圖21 有限元分析得到的管內壁等效塑性應變分布圖

對比圖20-a、20-b可知，解析模型預測的凹痕管道內壁的等效塑性應變的最大值均為0.160，位于凹痕中心，但應變分布范圍有所不同。圖21中有限元分析得到的管道內壁最大等效塑性應變?yōu)?.152，有限元結果與解析模型預測結果基本一致。因此可以認為，對于約束凹痕的總體應變評估，解析模型依然能很好地預測凹痕管道的最大應變值和峰值位置，且與文獻[21]所得到的結論一致。

4 結論

1）建立了含約束型和非約束型凹痕缺陷的有壓管道非線性有限元模型，分析了管道內外壁等效塑性應變沿管道軸向、環(huán)向的分布規(guī)律，結果表明對于內壓作用下初始加載深度為6%D（60.96 mm）的凹痕缺陷，整個凹痕管道的塑性化程度最大點位于凹痕最深處的管道內壁，與文獻所得到的結論一致，且約束凹痕的應變大于同等條件下的非約束凹痕。

2）通過對凹痕變形區(qū)域進行三維插值生成定義凹痕表面的位移函數(shù)，提出了以應變張量形式分析凹痕區(qū)域應變的解析模型，并利用該解析模型預測的應變分布與非線性有限元法預測的應變分布進行了對比，結果表明，解析模型可用來求解管道凹痕區(qū)域的應變，且能很好地預測凹痕管道的最大應變值和峰值位置。

3）通過考慮和不考慮剪切應變的等效塑性應變求解，發(fā)現(xiàn)剪切應變對凹痕區(qū)域等效塑性應變的峰值影響較小，因此在一般的等效塑性應變求解中剪切應變可以忽略，但剪切應變會對凹痕區(qū)域范圍的應變分布產生影響。所提出的解析模型可以高效便捷地獲得凹痕區(qū)域實際應變分布，為精確評價管道凹痕缺陷提供了科學依據，可用于改進現(xiàn)有的基于應變的凹痕評價體系。