基于中心線測量的山地管道力學分析軟件開發

蘭才富，趙 飛，趙 榮，吳承睿，郭家瑞，鄧寶信，王 琳

（1.西南管道南寧輸油氣分公司，廣西南寧 530200；2.西南石油大學，四川成都 610500）

0 引言

我國西南地區油氣管道沿線地質條件復雜、地災頻發，極易引起管道變形失效，因此需要對管道變形、失穩等狀態進行安全監測。然而，通過外部地質環境監測推測管道本體安全狀態的方法，沿線工作量大、難度高，且不能直接體現管道本體的真實變形；傳統的應力應變監測系統只能針對局部管段進行監測，全線應用的代價較大，而且多為事后監測，隱患監控滯后。因此，亟需開發對管道全線進行管道應力應變狀態計算、風險發展趨勢及安全評估的手段和工具[1-3]。

在此背景下，基于中心線測量結果，開發可對管道全線應變在空間上連續檢測的管道本體安全狀態評價方法：①利用單次測量結果，計算全線應力應變分布狀態；②利用多次測量結果，分析管道應力應變發展趨勢；③評價管道全線在不同工況下的可靠性，尋找出管道力學狀態的危險環節[4-6]。研究基于管道中心線檢測數據的管道應力應變算法，構建反映長輸管道全線應力應變狀態及其風險發展趨勢的管道安全評估方法，開發相應工具。研究結果可實現油氣管道的中心線、位移、應變應力的自主檢測和風險等級評估，為管道運行和維修維護決策提供適用、可靠、高效的技術工具，保障管道安全。

1 數學模型及求解方法

通過db4 小波的6 層軟閾值分解對角速度進行降噪處理，軟件的核心算法是管道力學分析的應力應變算法[7]。實現有限元分析標準化和規范化的載體為單元，其搭建出的復雜結構即為有限元分析的對象。有限元分析最主要的內容就是研究單元，首先給出單元的節點位移和節點力，然后基于單元節點位移與節點力的相互關系，獲得相應的剛度系數，進而得到單元的剛度方程，再將單元組裝為整體剛度方程，即可得到整體結構的基于節點位移的整體平衡方程。

基本分析流程為：節點編號和單元劃分→計算各單元剛度方程→組裝各單元剛度方程→處理邊界條件并求解→求支反力→求各單元其他力學量。

由于六面體單元的幾何形狀與管道結構相協調，因此與四面體相比，六面體構成環形單元的節點數目更少、網格質量更高，更便于在管道力學計算中實現對網格規模的控制。采用有限元法分析三維空間問題時，可將研究對象劃分為如圖1 所示的八節點六面體單元[8]。六面體單元模型由8 個節點組成，每個節點有3 個自由度，單元共計24 個自由度。

圖1 八節點六面體單元模型

其節點位移列陣qe與節點力列陣Pe分別表示為：

八節點單元的每個方向的位移場可以設定8 個待定系數，根據確定位移模式的基本原則（從低階到高階、唯一確定性），選取該單元的位移模式：

由節點條件確定待定系數（ai，bi，ci），i=0，1，2，……7。代入上式，可計算整理得到該單元的形狀函數矩陣：

得到該單元的形狀函數矩陣后，按照有限元分析的標準過程推導相應的幾何矩陣、剛度矩陣、節點等效載荷矩陣以及剛度方程。

由彈性力學平面問題的幾何方程，可得到單元應變場：

由彈性力學中平面問題的物理方程，可得到單元的應力表達，然后計算單元的勢能，就能得到單元的剛度矩陣與等效節點載荷矩陣：

將單元勢能對節點位移取一階極值，即可得到單元的剛度方程：

由六面體環形單元構造出管段的剛度矩陣K、彈性矩陣D、幾何矩陣B 后，進一步得管道總體的剛度方程：

管段總體位移向量后，由彈性力學空間問題的物理方程可以得到應力和應變場：

2 軟件設計理念及特色

為了保證軟件的實用性，詳細調研了國內外相關軟件（Abaqus、ANSYS）功能，明確軟件的設計目標為：吸收權威軟件的優勢功能，圍繞管道安全檢/監測管理的業務需求擴充（消弱）其部分功能，保證計算精度可靠的情況下，增強軟件的適用性；保證軟件操作簡便、直觀。為實現該目標，突破了常規力學分析軟件思路，引入新的軟件設計理念：采用模塊化編程，有機地吸取有限元軟件與信號分析軟件各自的優勢，實現計算軟件的精確性與大數據處理能力的有機結合。依據該理念，完成了西南管道基于中心線測量的山地管道力學分析軟件設計。

（1）具有較好的通用性，其數據預處理、中心線分析、特征識別、力學分析、安全評估模塊的功能可以單獨使用，滿足不同用戶的多方面需求。

（2）應用領域及用戶包括專業設計人員、專業技術人員、完整性管理人員及現場工人等用戶。

（3）使用邏輯與用戶需求相契合，滿足不同層次用戶需求，既可以通過導入數據進行快速計算得到結果，也可以用于對管道力學性能的詳細計算和分析。

（4）具有良好的人機界面，在自主開發的界面進行數據集成與處理，采用大量圖形和表格引導輸入數據和顯示計算結果，多頁面、積木式窗口設計，操作簡便，與各類常規軟件（Word、Excel、Access 等）之間可建立數據接口。

（5）采用高效的計算方法，計算效率要優于通用有限元軟件，并且操作界面簡潔、結果展示詳盡，具備一鍵報告生成功能。

（6）軟件界面有操作提示，并且對用戶的每一步操作都進行操作日志記錄。

3 軟件應用效果分析

軟件能否得以廣泛應用并發揮其效能，關鍵是中心線數據預處理、管道特征識別、應力應變算法的準確性。

3.1 數據預處理效果

小波變化作為一種有效的時頻分析方法，能夠在保留時間信息的前提下對不同頻率的信號進行分析[9-10]。“db”系小波在慣性導航的降噪上被廣泛的應用，通過“db4”小波6 層低頻重構得到的角速度信息與原始角速度信息對比如圖2 所示。

圖2 降噪前后角速度對比

3.2 特征識別方法使用效果

連續敷設管道采用彎頭來改變管道的走向，以滿足實際的運輸需求，在地質環境復雜的山區，彎頭的使用更為頻繁。為避免彎頭引起的姿態信息改變進而影響管道的應力應變計算，開發基于檢測數據的彎頭識別方法[11]。經過彎頭后，檢測器的姿態信息發生顯著變化，在角速度的信號中體檢為局部的明顯凸起，通過選取適當的閾值能夠有效的對直管和彎頭進行區分、通過曲率可以對彎頭的類型進行區分。

表1 為50 km 實際管道本軟件方法與第三方漏磁檢測的彎頭檢測結果對比，在50 km 管道，本方法共識別出722 處冷彎及342 處熱彎，共計1064 處彎頭，同管段第三方漏磁檢測結果為720 處冷彎、341 處熱彎；本方法的冷彎合計722 處、熱彎合計342 處，第三方檢測的冷彎合計720 處、熱彎合計341 處，誤差僅為0.65%。

表1 軟件識別彎頭及第三方檢測彎頭對比

3.3 應力應變算法驗證

采用本課題方法與Abaqus 計算結果進行對比驗證。管材彈性模量210 GPa，密度為7850 kg/m3，泊松比0.3，管徑為1016 mm，壁厚為16 mm，管長10 m。模擬集中載荷作用于水平管道中間頂部，表2 為不同工況下的集中載荷。

表2 不同工況下的集中載荷和最大應力

采用六面體環形單元計算得到的管道應力分布，與Abaqus 計算得到的管道應力分布結果吻合。

3.4 前50 km 安全評估結果

軟件的部分界面如圖3 所示，根據應力應變計算結果篩選得到的關注點共6 處，根據管道的應變均值及彎曲管道的長度，將彎曲管段分為局部彎曲和整段彎曲，長度劃分界限為2/3 選取管道長度。除了應力應變計算結果外，通過匹配內檢測數據還可以得到關注點的具體坐標信息，方便后續的實地考察及維護工作。

圖3 軟件界面

4 結論

基于MATLAB，開發了基于管道中心線檢測數據的管道應力應變評價軟件，能夠實現以下功能：

（1）對原始姿態信息進行降噪處理，以此提高應力應變計算結果的準確性。

（2）有效識別中心線數據中的彎頭數據，通過中心線數據能夠識別彎頭的類型，曲率半徑及彎頭的走向。

（3）采用六面體法沿線計算管道的應力應變。

（4）通過應力應變計算結果、彎曲變形管段長度對發生變形的危險管道進行識別及分類。

通過與原始姿態信號、第三方檢測結果、Abaqus仿真計算的對比，分別驗證了該軟件數據預處理模塊、特征識別模塊、應力應變計算模塊的正確性。該軟件可通過中心線檢測結果，快速給出全線的管道應力應變水平，及時發現管道中有應力應變集中的管段，為長輸管道的安全評估提供有效的方法。