基于邊界元法的站場陰極保護數值模擬軟件設計與開發

廖柯熹,馬曉磊,張中放,周衛軍,張 瑤,任馳遠(. 西南石油大學 石油與天然氣工程學院，成都 60500； . 中國石油天然氣集團公司 塔里木油田分公司油氣運銷部，庫爾勒 84000； . 西南石油大學 理學院，成都 60500)

基于邊界元法的站場陰極保護數值模擬軟件設計與開發

廖柯熹1,馬曉磊1,張中放2,周衛軍2,張 瑤2,任馳遠3
(1. 西南石油大學 石油與天然氣工程學院，成都 610500； 2. 中國石油天然氣集團公司 塔里木油田分公司油氣運銷部，庫爾勒 841000； 3. 西南石油大學 理學院，成都 610500)

邊界元法在陰極保護電位分布計算方面有著明顯的優勢，采用邊界元法對數值模型邊界進行離散，計算節點少，計算速率快。站場陰極保護數值模擬軟件基于邊界元法對站場陰極保護電位分布進行計算，得到可視化結果，通過該結果可以發現屏蔽和干擾位置，對輔助陽極位置進行優化，從而解決站場陰極保護系統腐蝕區域不能預測的困難，使站場陰極保護系統工程更加合理有效。

邊界元法；站場；陰極保護；電位分布模擬

20世紀60年代國外開始進行區域性陰極保護的研究和應用，我國從20世紀70年代末、80年代初開始在油田和部分輸油站嘗試采用區域性陰極保護技術，到90年代中期該技術已相對成熟并開始推廣應用。2001年11月，鄯烏輸氣管道，鄯善首站成為第一個實施區域性陰極保護的壓氣站[1]。

目前，現場陰極保護施工過程中普遍采用試湊法，保護效果不太理想，往往需要進行后期整改。隨著數值模擬技術在陰極保護領域的深入應用，使得復雜環境中陰極保護效果的預測和低成本分析研究成為可能。與傳統方法相比，數值模擬技術在管道和站場陰極保護各種影響因素的干擾趨勢和規律研究方面具有無法比擬的優勢，同時可用于陽極位置優化、陰極保護效果的預測和全面評價[2]。

目前，可以進行區域陰極保護系統電位分布計算的軟件主要是基于有限元法和邊界元法進行計算，基于有限元法進行計算的軟件有FLUENT、COMSOL等，網格劃分復雜，計算節點多，不適用于現場情況。而BEASY是基于邊界元法進行計算的軟件，在進行網格劃分時，只需對邊界進行劃分，計算量小，可以節省計算時間。

由于陰極保護計算通常只關心材料外表面和地表的參數分布情況，相較于有限元法，邊界元法具有只對邊界進行離散化的特點，更適用于計算陰極保護等與腐蝕相關的問題。為此，本工作設計開發了基于邊界元法的陰極保護電位分布計算軟件——站場陰極保護數值模擬軟件(Numerical Simulation Software Of Cathodic Protection，簡稱NSSCP)，以期可靠分析與評價站場陰極保護系統的有效性。

1 陰極保護電位數值計算模型

埋地管道陰極保護電位控制方程和儲罐底板陰極保護電位控制方程的推導過程類似，在此僅考慮儲罐底板陰極保護的電位控制方程，可以用Poisson方程表示為

2φ

式中：δ(X-Xe)表示以陽極點X(x,y,z)為中心的狄拉克函數，其性質見式(2)。

式中：φ是施加陰極保護后的極化電位值；x，y和z分別是計算點與儲罐底板中心的距離；Ie為陽極點Xe(x,y,z)輸出的電流(即點源強度)。當所研究的恒定電場中存在恒定的電源，電位滿足poisson方程，在給定邊界條件時，可通過求解該方程得到電位的分布。

公式(1)僅用來描述儲罐底板陰極保護的電位分布，但是如果想要求得具體解，就需要附加邊界條件。數值模型的邊界為Ω(研究區域)和總邊界Γ(Γ表面包圍3 Ω),總邊界Γ包括無限大地面表面邊界Γd、儲罐底板外表面邊界Γc、虛設半無限域土壤球冠邊界?！蕖?/p>

邊界條件簡化：地表面邊界Γd與大氣相接觸，外加陰極保護電流只能在土壤介質中流動，不能流入大氣，所以把地表面當作絕緣面；假設無窮遠處的電位為零，忽略陽極對半無限大土壤表面的影響，所以在無窮遠處表面的電流密度也為零；陰極的邊界條件采用描述電流密度J和電位φ之間函數關系的極化函數來定義

J=f(φ

式中：φeq是儲罐底板在土壤環境中的自腐蝕電位，f(φ-φeq)為陰極極化函數。

綜上所述得到數學模型及邊值問題如下

2 邊界積分方程推導與處理

邊界元方法的求解是用格林公式對微分方程積分后再離散處理，最終通過求解關于電位和電流密度的線性方程組系統?？杀硎緸?/p>

式中：H、G分別為通過邊界元方法定義的系數矩陣；E、I分別為節點電位矩陣和電流密度矢量，其矩陣中每一項分別為待求的電位和電流密度[3]。

2.1 邊界積分方程推導

以P′為圓心，以微小ε為半徑，在邊界上建立微小半球模型，設Q′點位于邊界上。通過積分方程式把域內的點用P、Q表示，邊界上的點用P′、Q′的位勢值φ以及通量值q聯系起來，對區域Ω內任意點均成立。在陰極保護電位分布模型中，邊界積分方程式為

C(P′)φ(P′)

所研究的模型屬于半無限域，很難在虛設半無限域土壤球冠邊界Γ∞和無限大地面Γd表面邊界上劃分網格。因此根據模型已知的邊界條件，對模型的邊值問題進行簡化。在邊界Γ∞處，滿足函數值和函數的法向導數都為零，即：

代入后可以發現，不需要在邊界?！奚蟿澐謫卧陀嬎?。在邊界Γd上：

代入后不能完全消除該邊界的積分項，因此需要對邊界Γd進行單元劃分。由于Γd是邊界無限大地面表面，不能確定應該取多大的邊界來劃分單元，同時劃分單元后也會增大計算量。通過鏡像法可在Γd邊界，尋找一個特定的基本解φ*，使邊界Γd無需劃分。

根據Poisson方程的基本解得到的模型為

φ*(b)=φ*(A,B)+φ*(A′,B′)=

根據模型基本解，邊界積分方程可簡化為

選用模型的基本解時，無需對邊界Γd劃分單元，地表面的邊界條件自動滿足，所以只需對陰極的邊界進行劃分。

2.2 邊界離散

邊界元法將邊界上的積分方程式離散化，并歸結為代數方程組求得近似解。通過插值函數得到各個邊界單元上的函數值和函數的法向導數值。

按單元求和的離散形式改寫為按節點求和的邊界積分方程離散格式，即

由于j點也是邊界上的點，所以必然與某I節點為同一點，所以合并同類項后，公式簡化為

式(14)中:

對于m個節點，可以得到聯立的一次方程組為

式中:H與G表示系數矩陣，Φ表示邊界單元節點處的函數值的列向量，Q表示邊界單元節點處函數法向導數值的列向量，按照已給定的邊界條件，由式(16)可以計算出邊界上節點的未知函數值與其法向導數值。

2.3 輔助陽極處理

根據δ函數的性質可知，域內存在一個集中源點X1時

式中:φ*(P′,A1)是對應A1的基本解，得

當域內存在多個集中源點An(n=1LN)時，得

φ*(P′,An)b(Q)dΩ(Q)=

積分項Bj不會產生未知量。

2.4 非線性邊界條件處理

在數值模擬過程中，邊界條件的確定及其處理方式直接影響計算結果的準確性，軟件總共有4種邊界條件，包括極化曲線邊界、電位邊界、電流邊界、阻抗邊界，其中最主要的邊界條件是非線性極化曲線邊界。

由于陰極和陽極的電位和電流密度的關系并非簡單的線性關系，而是非線性關系，在模型處理過程中必須對其進行處理。儲罐底板表面受到多種不同因素的影響，極化函數不能只用單獨的環境或者材料的影響函數來表達，為得到極化函數定義的陰極邊界條件J=f(φ-φeq)，通過試驗測得極化曲線作為邊界條件是非常有必要的。

通過對測得的極化曲線進行擬合，可以確定電流密度與極化電位之間的定量關系。主流的擬合方法有Newton-Raphson迭代法和分段線性擬合法。

由于被保護金屬表面受到多種不同因素的影響，所以極化函數不能只用單獨的環境或者材料的影響函數來表達，為得到極化函數定義的陰極邊界條件，綜合考慮現場實際狀況，通過試驗測得被保護金屬在土壤中的電位分布，為了取得良好的收斂性和計算精度，該軟件采用分段擬線性化法處理非線性邊界條件。擬合得到電流密度J與電位φ的之間的線性關系

J=aφ

推導得

將式(22)代入下列矩陣方程

移項整理得

矩陣最終化為AX=F的形式，輔以迭代的方法，可以求解得到電位值大小。

3 軟件開發

基于邊界元算法開發了“站場陰極保護數值模擬軟件(NSSCP)”，該軟件可以實現三維站場陰極保護系統的模擬，獲得可視化電位分布規律，從而解決現有站場和新改擴建站場內埋地金屬構筑物密集所造成的屏蔽和干擾問題，并防止站場內陰極保護電流對未納入陰極保護管道的干擾影響。相較于有限元法計算,該軟件的優點就是簡單快捷。

3.1 軟件開發流程和計算流程

在軟件的設計開發過程中，使用Access 2007數據庫作為后臺數據庫，并基于Visual Studio C# 2010設計了軟件用戶界面。軟件共劃分出陰極保護系統參數存儲、邊界條件處理、電位計算及數據導出四個主模塊。圖1和圖2分別是軟件開發的流程圖和軟件計算流程圖。

圖1 軟件開發流程Fig. 1 Software development process

圖2 軟件計算流程Fig. 2 Software calculation process

3.2 陰極保護系統參數存儲和建模

用戶可以在“工程參數”界面輸入工程名稱、站場名稱、土壤電阻、站場簡介、模擬目的、精度選擇(高精度、一般精度、較低精度)、精度參數(最大迭代次數和收斂誤差閾值)等基礎數據。土壤電阻率會隨地域、季節和天氣發生變化，在設置土壤電阻率時，應考慮這些因素的影響。在“實體對象”界面中可以進行腐蝕實體對象建模，儲罐、陽極、埋地管道需要設置保護層系數、直徑、埋深、位置坐標、邊界條件等參數。

3.3 電位分布模擬計算

采用對于求解區域進行自動剖分網格的方法來劃分單元，將形成的邊界單元的坐標保存在相應的輸入文件中。計算系數矩陣，調用邊界條件，求解線性方程組，完成迭代。

軟件在使用時最重要的步驟就是進行參數收集與導入，在軟件的主界面應先進行陰極極化曲線邊界參數的導入；其次是土壤電阻率的導入，需要考慮不同季節導致的土壤電阻率的變化，根據實測結果進行導入；最后導入的是實體腐蝕對象的參數(包括地理位置、幾何尺寸等)，通過一系列的計算便可以得到電位分布。

干擾分析菜單欄可以實現對計算結果查看和分析，得到整個站場節點電位分布數據列表,見圖3。

干擾問題是站場陰極保護系統實施的一個難題。

圖3 數值計算結果查看界面Fig. 3 View interface of numerical results

對油氣輸送站場進行陰極保護，由于是將整個站區內的所有地下金屬結構全部納入保護系統，因此一般不會產生內部干擾，但是由于站場陰極保護電流通常遠大于干線保護電流，故常常對管道干線等外部結構及其陰極保護系統造成干擾。一般系統電流輸出越大，干擾也會越嚴重。通過該軟件干擾分析產生的干擾數據可以有效判斷出是否發生發生干擾，提高了站場陰極保護系統的可靠性。

4 實例驗證

某站場使用一套陰極保護系統對四個儲罐進行保護，6組深井陽極埋設45 m分別位于儲罐的兩側每組陽極由8支研究塊組成,共8 m,直徑為25 cm，儲罐直徑為60 m，恒電位儀的輸出電位為8 V。對儲罐周邊測試樁進行檢測，測得電位均處于-0.85～-1.2 V(相對銅/硫酸銅參比電極)，達到保護狀態，保護效果良好。

通過采用COMSOL軟件模擬結果進行比對分析，模擬時，所取土壤電阻率為50 Ω·m，砂墊層200 Ω·m，瀝青砂層1 000 Ω·m。采用COMSOL軟件模擬得到的儲罐保護電位分布云圖，通過模擬計算，并將COMSOL軟件模擬結果和實測數據進行統計,見表1。

表1 計算結果比較Tab. 1 Comparison of calculated results

COMSOL基于有限元方法進行計算，計算時間較長，一般認為其結果較為精確，但是存在外邊界截斷誤差；而NSSCP是基于邊界元法進行開發的，不存在邊界截斷誤差。由表1可見，NSSCP的計算精度接近COMSOL的，證明該軟件計算結果可靠，可以作為站場陰極保護系統的檢測維護軟件。

5 結束語

(1) NSSCP軟件在計算結果精度方面接近COMSOL軟件，并且計算過程中收斂速率快，計算時間短，提高了使用效率。

(2) NSSCP軟件具有較高的實用性，實際應用效果較好，能夠用于指導站場陰極保護系統的設計計算和科學維護。

[1] 陳洪源,范志剛,劉玲莉,等. 區域性陰極保護技術在輸氣站場中的應用[J]. 油氣儲運,2005,24(5):41-44.

[2] 張豐,陳洪源,李國棟,等. 數值模擬在管道和站場陰極保護中的應用[J]. 油氣儲運,2011,30(3):208-212.

[3] 陳石,秦朝葵. 邊界元法介紹及其在陰極保護領域應用[J]. 上海煤氣,2015(3):7-10.

“2017材料物理測試新技術研討會”征 文 通 知

由上海材料研究所等單位聯合主辦，《理化檢驗-物理分冊》編輯部、《機械工程材料》編輯部共同承辦的“物理測試新技術研討會”已連續成功舉辦兩屆，受到了國內廣大材料物理測試專家、學者及從業人員的熱烈歡迎和積極參與。在此基礎上，“2017材料物理測試新技術研討會”將繼續圍繞“新的測試技術與方法”、“新的測試儀器和設備”、“新的試驗標準與規范”等方面展開討論，旨在為廣大物理測試工作者以及設備廠商提供一個溝通、交流和展示的平臺，以總結經驗，促進創新，進而推動我國材料物理測試技術的發展與進步。本次會議將于2017年10月底在上海舉行，由上海材料研究所、中國機械工程學會理化檢驗分會、中國機械工程學會材料分會聯合主辦，《理化檢驗-物理分冊》編輯部、《機械工程材料》編輯部和上海市工程材料應用與評價重點實驗室共同承辦。會議將邀請業內權威專家作專題報告，同時熱忱歡迎行業專家、學者及從業人員踴躍投稿，積極參會。所征論文將以《理化檢驗-物理分冊》正刊專欄或增刊形式發表，評選出的優秀論文將在會議期間進行口頭交流。

1 征文內容

(1) 材料物理測試的新技術、新方法、新標準和新儀器設備的介紹及應用；

(2) 物理測試技術在新材料開發、研究中的應用；

(3) 物理測試技術在產品質量控制過程中的應用；

(4) 物理測試技術在材料表面處理工藝、材料冷熱加工工藝、材料改性技術以及先進制造技術方面的應用；

(5) 物理測試技術在產品在役安全監控、壽命估算、可靠性評估等方面的應用；

(6) 機械產品及其構件失效分析的最新思路、技術、方法及典型的失效分析案例介紹；

(7) 試驗設備的升級改造技術；

(8) 試驗室的管理及在全面質量管理中的成功經驗；

(9) 其他與材料物理測試技術相關的論文。

2 征文要求

(1) 未公開發表過的原創性論文。

(2) 來稿要求論點明確、數據可靠、邏輯嚴密、文字精煉。每篇論文必須包括題目、作者姓名、作者單位、單位所在地及郵政編碼、摘要和關鍵詞、正文、參考文獻和第一作者及通信作者簡介(包括姓名、性別、職稱、出生年月、所獲學位、目前主要從事的工作和研究方向)，在文稿的首頁地腳處注明論文屬何項目、何基金(編號)資助，沒有的不注明。論文題目、作者姓名、作者單位、單位所在地及郵政編碼、摘要和關鍵詞應中英文對照。

(3) 論文摘要盡量寫成報道性文摘，包括目的、方法、結果、結論4方面內容(200字左右)，應具有獨立性與自含性，關鍵詞選擇貼近文義的規范性單詞或組合詞(3～8個)。

(4) 文稿篇幅(含圖表)一般不超過8 000字。文中量和單位的使用請參照中華人民共和國法定計量單位最新標準。外文字符必須分清大、小寫，正、斜體，黑、白體，上下角標應區別明顯。

(5) 文中的圖、表應有自明性，圖像要清晰，層次要分明；所有圖和表均應標明中英文圖題和表題。

(6) 參考文獻的著錄格式采用順序編碼制，請按文中出現的先后順序編號。所引文獻必須是作者直接閱讀參考過的、最主要的、公開出版文獻。未公開發表的、且很有必要引用的，請采用腳注方式標明，參考文獻不少于5條。

3 征文截止時間

2017年8月31日。

4 論文提交

論文投稿至E-mail: pt@mat-test.com，請注明稿件為“物理測試會議論文”。

5 聯系方式

地址： 上海市邯鄲路99號，《理化檢驗-物理分冊》編輯部、《機械工程材料》編輯部。

聯系人：李玲；金靜靜。

電話：021-65559079；021-65556775-361。E-mail：pt@mat-test.com。

(會議的具體時間與地點，及會務費用等詳情將在第二輪通知中發出。)

Design and Development of Numerical Simulation Software for Cathodic Protection of Stations Based on Boundary Element Method

LIAO Kexi1, MA Xiaolei1, ZHANG Zhongfang2, ZHOU Weijun2, ZHANG Yao2, REN Chiyuan3
(1. College of Petroleum & Natural Gas Engineering, Southwest Petroleum University, Chengdu 610500, China;2. Tarim Oilfield Company Oil and Gas Distribution Unit, China National Petroleum Corporation, Korla 841000, China;3. College of Science, Southwest Petroleum University, Chengdu 610500, China)

Boundary element method in distributed computing of cathodic protection potential has a distinct advantage. When the boundary element method to scatter a numerical model was used, the calculation nodes were less and calculation speed was fast. The software based on the boundary element method calculated the cathodic protection potential distribution in cathodic station, and the results were visual, shielding and interference position could be found through the results, auxiliary anode position could be optimized in order to resolve the difficulty that the station cathodic corrosion protection system area could not be predicted, so the station cathodic protection system engineering could be more rational and effective.

boundary element method; station; cathodic protection; potential distribution simulation

2016-01-15

廖柯熹(1970-)，教授，博士，從事油氣儲運系統完整性管理與安全評價工作,liaokxswpi@163.com

10.11973/fsyfh-201707015

TE88

A

1005-748X(2017)07-0551-06