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

基于邊界元法的站場陰極保護(hù)數(shù)值模擬軟件設(shè)計與開發(fā)

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

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

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

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

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

目前，可以進(jìn)行區(qū)域陰極保護(hù)系統(tǒng)電位分布計算的軟件主要是基于有限元法和邊界元法進(jìn)行計算，基于有限元法進(jìn)行計算的軟件有FLUENT、COMSOL等，網(wǎng)格劃分復(fù)雜，計算節(jié)點多，不適用于現(xiàn)場情況。而BEASY是基于邊界元法進(jìn)行計算的軟件，在進(jìn)行網(wǎng)格劃分時，只需對邊界進(jìn)行劃分，計算量小，可以節(jié)省計算時間。

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

1 陰極保護(hù)電位數(shù)值計算模型

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

2φ

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

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

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

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

J=f(φ

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

綜上所述得到數(shù)學(xué)模型及邊值問題如下

2 邊界積分方程推導(dǎo)與處理

邊界元方法的求解是用格林公式對微分方程積分后再離散處理，最終通過求解關(guān)于電位和電流密度的線性方程組系統(tǒng)。可表示為

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

2.1 邊界積分方程推導(dǎo)

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

C(P′)φ(P′)

所研究的模型屬于半無限域，很難在虛設(shè)半無限域土壤球冠邊界?！藓蜔o限大地面Γd表面邊界上劃分網(wǎng)格。因此根據(jù)模型已知的邊界條件，對模型的邊值問題進(jìn)行簡化。在邊界Γ∞處，滿足函數(shù)值和函數(shù)的法向?qū)?shù)都為零，即：

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

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

根據(jù)Poisson方程的基本解得到的模型為

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

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

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

2.2 邊界離散

邊界元法將邊界上的積分方程式離散化，并歸結(jié)為代數(shù)方程組求得近似解。通過插值函數(shù)得到各個邊界單元上的函數(shù)值和函數(shù)的法向?qū)?shù)值。

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

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

式(14)中:

對于m個節(jié)點，可以得到聯(lián)立的一次方程組為

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

2.3 輔助陽極處理

根據(jù)δ函數(shù)的性質(zhì)可知，域內(nèi)存在一個集中源點X1時

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

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

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

積分項Bj不會產(chǎn)生未知量。

2.4 非線性邊界條件處理

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

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

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

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

J=aφ

推導(dǎo)得

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

移項整理得

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

3 軟件開發(fā)

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

3.1 軟件開發(fā)流程和計算流程

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

圖1 軟件開發(fā)流程Fig. 1 Software development process

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

3.2 陰極保護(hù)系統(tǒng)參數(shù)存儲和建模

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

3.3 電位分布模擬計算

采用對于求解區(qū)域進(jìn)行自動剖分網(wǎng)格的方法來劃分單元，將形成的邊界單元的坐標(biāo)保存在相應(yīng)的輸入文件中。計算系數(shù)矩陣，調(diào)用邊界條件，求解線性方程組，完成迭代。

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

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

干擾問題是站場陰極保護(hù)系統(tǒng)實施的一個難題。

圖3 數(shù)值計算結(jié)果查看界面Fig. 3 View interface of numerical results

對油氣輸送站場進(jìn)行陰極保護(hù)，由于是將整個站區(qū)內(nèi)的所有地下金屬結(jié)構(gòu)全部納入保護(hù)系統(tǒng)，因此一般不會產(chǎn)生內(nèi)部干擾，但是由于站場陰極保護(hù)電流通常遠(yuǎn)大于干線保護(hù)電流，故常常對管道干線等外部結(jié)構(gòu)及其陰極保護(hù)系統(tǒng)造成干擾。一般系統(tǒng)電流輸出越大，干擾也會越嚴(yán)重。通過該軟件干擾分析產(chǎn)生的干擾數(shù)據(jù)可以有效判斷出是否發(fā)生發(fā)生干擾，提高了站場陰極保護(hù)系統(tǒng)的可靠性。

4 實例驗證

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

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

表1 計算結(jié)果比較Tab. 1 Comparison of calculated results

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

5 結(jié)束語

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

(2) NSSCP軟件具有較高的實用性，實際應(yīng)用效果較好，能夠用于指導(dǎo)站場陰極保護(hù)系統(tǒng)的設(shè)計計算和科學(xué)維護(hù)。

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

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

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

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

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

1 征文內(nèi)容

(1) 材料物理測試的新技術(shù)、新方法、新標(biāo)準(zhǔn)和新儀器設(shè)備的介紹及應(yīng)用；

(2) 物理測試技術(shù)在新材料開發(fā)、研究中的應(yīng)用；

(3) 物理測試技術(shù)在產(chǎn)品質(zhì)量控制過程中的應(yīng)用；

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

(5) 物理測試技術(shù)在產(chǎn)品在役安全監(jiān)控、壽命估算、可靠性評估等方面的應(yīng)用；

(6) 機(jī)械產(chǎn)品及其構(gòu)件失效分析的最新思路、技術(shù)、方法及典型的失效分析案例介紹；

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

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

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

2 征文要求

(1) 未公開發(fā)表過的原創(chuàng)性論文。

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

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

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

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

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

3 征文截止時間

2017年8月31日。

4 論文提交

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

5 聯(lián)系方式

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

聯(lián)系人：李玲；金靜靜。

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

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

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-)，教授，博士，從事油氣儲運系統(tǒng)完整性管理與安全評價工作,liaokxswpi@163.com

10.11973/fsyfh-201707015

TE88

A

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