羅 威,袁銳波,陳有錦,朱 正,唐世旭,胡啟明
(1.昆明理工大學(xué) 機(jī)電工程學(xué)院, 昆明 650504;2.中石化勝利油田分公司海洋采油廠, 山東 東營(yíng) 257237)
隨著社會(huì)經(jīng)濟(jì)發(fā)展,人們對(duì)物質(zhì)需求的不斷提升,世界各國(guó)對(duì)石油資源的需求日益增長(zhǎng),其中,我國(guó)石油年消費(fèi)量由2001年的2.28億噸上升至2020年的6.54億噸。隨著石油消費(fèi)日益增加,管道運(yùn)輸作為石油運(yùn)輸?shù)闹饕绞街?在世界范圍受到廣泛應(yīng)用,據(jù)統(tǒng)計(jì),我國(guó)2020年油氣管道總里程已達(dá)1.44×105km,預(yù)計(jì)在2025年,我國(guó)油氣管網(wǎng)總長(zhǎng)將達(dá)到2.4×105km。對(duì)于如此長(zhǎng)的管道里程,如何保證它們安全穩(wěn)定運(yùn)行非常值得研究。
隨著石油管道使用時(shí)間增加,管道的腐蝕、磨損和老化情況日益嚴(yán)重,容易造成泄漏,安全隱患極大。如果石油管道發(fā)生泄漏,不僅管道輸送效率的降低,更會(huì)對(duì)周圍環(huán)境造成污染,引發(fā)一系列事故,威脅到人民群眾的生命安全[1]。因此,當(dāng)石油管道發(fā)生泄漏時(shí),如果能及時(shí)發(fā)現(xiàn)泄漏并對(duì)泄漏位置進(jìn)行定位,將大大降低石油管道泄漏的危害程度。
基于模型檢測(cè)的管道檢測(cè)方法主要是根據(jù)管道的流體動(dòng)力學(xué)建立模型或其他方法建立整個(gè)管道模型,通過(guò)已知的邊界條件,對(duì)管道沿線的壓力、流量和溫度等流體參數(shù)進(jìn)行預(yù)估計(jì),將預(yù)估值與實(shí)際值進(jìn)行比較,以此判斷管道是否發(fā)生泄漏并實(shí)現(xiàn)定位功能[2-3]。由于泄漏位置不確定,該方法需要在管道不同位置選擇泄漏點(diǎn)進(jìn)行計(jì)算,依據(jù)龐大的計(jì)算來(lái)實(shí)現(xiàn)泄漏定位。為了在實(shí)現(xiàn)定位的同時(shí)減少計(jì)算量,將優(yōu)化算法與基于模型檢測(cè)方法相結(jié)合是很多學(xué)者的選擇方向,如Pudar等[4]選擇Levenberg-Marquardt算法求解管道泄漏模型,但該算法對(duì)初值極為敏感,容易陷入局部最優(yōu)解,同時(shí)該算法在迭代過(guò)程需要耗費(fèi)大量時(shí)間計(jì)算海森矩陣;Kapelan等[5]總結(jié)了前人的不足,將遺傳算法與Levenberg-Marquardt算法相結(jié)合,先是利用遺傳算法全局搜索,然后用Levenberg-Marquardt算法局部精準(zhǔn)搜索,但是這種組合方式會(huì)降低計(jì)算效率;陳特歡等[6]利用粒子群算法求解管道泄漏模型,并對(duì)參數(shù)敏感性進(jìn)行了分析,但是該算法需要迭代許多次才能收斂到較好的狀態(tài);華東陽(yáng)等[7]提出了一種基于遺傳算法的求解泄漏管道模型方法,該方法具備良好的準(zhǔn)確性和適應(yīng)性,但只考慮壓力的求解模型不夠準(zhǔn)確。
針對(duì)前面提到的國(guó)內(nèi)外各學(xué)者在基于模型法的管道檢測(cè)上出現(xiàn)的收斂速度慢,管道泄漏模型不夠準(zhǔn)確等問題,本文建立一個(gè)更為準(zhǔn)確的管道泄漏模型,并利用鯨魚優(yōu)化算法對(duì)其進(jìn)行求解,還對(duì)鯨魚優(yōu)化算法、粒子群算法、遺傳算法進(jìn)行多方面對(duì)比,以體現(xiàn)鯨魚優(yōu)化算法的優(yōu)勢(shì)。
根據(jù)動(dòng)量守恒、質(zhì)量守恒和能量守恒三大定律,可以推算出管道流體流動(dòng)過(guò)程中動(dòng)量方程、連續(xù)性方程和能量方程,將這3個(gè)方程組成偏微分方程組,可基本概括管道流體流動(dòng)的過(guò)程,該偏微分方程組就是流體管道的動(dòng)力學(xué)模型[8],可以寫成如下形式:
(1)
式中:F為自變量的向量,可以表示為F=[P,V,T]T,其中P為管內(nèi)壓力,Pa;V為流體流速,m/s;T為流體溫度,K;N為系數(shù)矩陣;M是與系數(shù)無(wú)關(guān)項(xiàng)的列向量;t、x分別表示時(shí)間變量和管道位置變量。
管道泄漏屬于瞬變狀態(tài),泄漏時(shí)管內(nèi)流體的切應(yīng)力τ0有別于定常流,因此,需要對(duì)切應(yīng)力τ0進(jìn)行修正,考慮本地加速度和對(duì)流項(xiàng)的影響[9-11],其修正結(jié)果為:
(2)
式中:f為摩擦因數(shù),取絕對(duì)值是為了保證切應(yīng)力方向始終與速度方向相反;d為管道內(nèi)徑,m;kt、kx分別是本地加速度系數(shù)和對(duì)流項(xiàng)系數(shù);sign(V)為符號(hào)項(xiàng),當(dāng)速度大于等于零時(shí),符號(hào)項(xiàng)為正,反之符號(hào)項(xiàng)為負(fù);a為流體波速,m/s。
因?yàn)槭凸艿涝阡佋O(shè)時(shí)一般選擇平直鋪設(shè),減少高程變化帶來(lái)的影響,因此,本文只探討沒有高程變化的平直管道。對(duì)于沒有高程變化的石油管道,根據(jù)修正后的切應(yīng)力方程,可以將管道模型中的N和M寫成如下形式:
(4)
式中:ρ為流體密度,kg/m3;Cv為定壓比熱容,J/(kg·K);D為管道外徑,m;K為管壁傳熱系數(shù),W/(m·K);T0為環(huán)境溫度,K。
上述的雙曲型偏微分方程組是不存在通解的,但是可以通過(guò)特征線法對(duì)其求解。根據(jù)方程組的3條特征線將其轉(zhuǎn)化為特殊的全微分方程,然后將這些全微分方程分別沿著它們的特征線進(jìn)行積分,得到更好處理的有限差分方程。首先根據(jù)線性代數(shù)理論,可以將管道模型寫成如下的全微分方程:
(5)
式中:
(6)
(7)
圖1是特征線法求解管道模型的示意圖,圖中:O、Q、G、H表示不同空間和時(shí)間位置節(jié)點(diǎn)的流動(dòng)參數(shù)壓力P、速度V和溫度T的關(guān)系,tj時(shí)刻S、R、V點(diǎn)所在位置的流動(dòng)參數(shù)可以通過(guò)O、Q、G的流動(dòng)參數(shù)獲得,然后將全微分方程積分,可以獲得tj+1時(shí)刻H位置的流動(dòng)參數(shù),這樣就可以沿著時(shí)間方向依次計(jì)算各個(gè)節(jié)點(diǎn)的流動(dòng)參數(shù)。
圖1 特征線法求解示意圖
將全微分方程分別沿著它們的特征線C+、C、C-進(jìn)行積分。積分結(jié)果如下:
PH=(PS-PV)+
(8)
(9)
TH=TR-
(10)
由此,在已知管道前一時(shí)刻的流動(dòng)參數(shù)下,可以通過(guò)方程(8)—(10)計(jì)算這一時(shí)刻的流動(dòng)參數(shù),因此,只需要知道管道初始的壓力、速度和溫度,就能對(duì)整個(gè)流體管道進(jìn)行求解。
對(duì)于流體管道模型的邊界條件來(lái)說(shuō),其可選范圍非常廣,只需要知道壓力、速度和溫度在輸入端或輸出端的值就行,所以邊界條件的組合一共有8種。本文根據(jù)實(shí)際情況選擇的邊界條件為:已知輸入端的壓力、溫度和輸出端的速度。
對(duì)于泄漏管道而言,其管道內(nèi)流體的流動(dòng)特性依然可以用上一小節(jié)的管道模型描述,但是在泄漏位置處,流體的流動(dòng)特性會(huì)發(fā)生改變。因此,對(duì)于泄漏管道模型的構(gòu)建,只需要在原來(lái)的模型上添加一個(gè)特殊的泄漏邊界條件即可[12-13]。
在如圖2所示的管道泄漏示意圖中,流體從左向右流動(dòng),并且在泄漏口周圍設(shè)置了一個(gè)小區(qū)域,稱為“泄漏區(qū)域”,針對(duì)這個(gè)區(qū)域設(shè)定如下的特殊邊界條件。假設(shè)泄漏孔很小,當(dāng)管道發(fā)生泄漏時(shí),可以認(rèn)為泄漏孔所在的截面上壓力變化極小,即各個(gè)方向上的壓力都是相同的,所以能得到下式:
圖2 管道泄漏示意圖
PH,i-=Pk=PH,i+
(11)
另外泄漏孔的泄漏量可以根據(jù)小孔方程求得:
(12)
式中:Ck、Ak分別為管道泄漏孔的流通系數(shù)和流通面積。根據(jù)流體的連續(xù)性還能得到如下方程:
Qk=AVH,i--AVH,i+
(13)
聯(lián)立式(8)—(13)求解得:
(14)
上式中:
根據(jù)前面所述的管道模型求解方式,使用Matlab編程語(yǔ)言對(duì)單管石油管道的瞬態(tài)模型求解。為了驗(yàn)證切應(yīng)力修正后管道模型的準(zhǔn)確性,還使用Pipeline Studio軟件對(duì)同樣的管道進(jìn)行仿真分析,然后將仿真數(shù)據(jù)和Matlab數(shù)據(jù)對(duì)比分析。
本次仿真數(shù)據(jù)來(lái)源于國(guó)內(nèi)某段油氣管道,該管道全長(zhǎng)84.7 km,由多個(gè)管段組成,管道材質(zhì)為L(zhǎng)450直縫鋼管,管道內(nèi)徑為914 mm,粗糙度為25.4,壁厚為12.7 mm,設(shè)計(jì)壓力為8.5 MPa,當(dāng)輸入端壓力為7 MPa,輸入端溫度為56 ℃,輸出端流量為3 500 m3/h,泄漏點(diǎn)在距輸入端65 km處時(shí),自編程序和仿真軟件計(jì)算的管道沿程壓力、速度如圖3和圖4所示。
圖3 自編程序和仿真軟件的管道沿程壓力曲線
圖4 自編程序和仿真軟件的管道沿程壓力曲線
自編程序和Pipeline Studio計(jì)算的沿程壓力最大相對(duì)誤差為1.02%,計(jì)算的沿程流速最大相對(duì)誤差為3.84%,均在工程誤差接受范圍內(nèi)。誤差的主要來(lái)源有2部分,一部分是由于分段計(jì)算時(shí)每段保留的有效數(shù)字不同,誤差會(huì)逐段累計(jì);還有一部分誤差來(lái)自于泄漏處邊界條件計(jì)算公式的差異。由此驗(yàn)證切應(yīng)力修正后管道模型具有較高的精度。
對(duì)整個(gè)管道而言,發(fā)生泄漏時(shí)會(huì)引發(fā)一系列影響,而這些影響將是判斷是否發(fā)生泄漏和定位的依據(jù)。對(duì)石油管道來(lái)說(shuō),泄漏后管道內(nèi)的工況將發(fā)生以下變化[14]:① 壓力下降:石油泄漏會(huì)導(dǎo)致管道內(nèi)部的壓力下降,因?yàn)楣艿劳ㄟ^(guò)泄漏孔與外界大氣壓相連。② 流速變化:泄漏會(huì)導(dǎo)致管道內(nèi)石油的流速發(fā)生變化。在泄漏口處,流速會(huì)很大,隨著距離泄漏口越來(lái)越遠(yuǎn),流速會(huì)逐漸降低。③ 溫度變化:當(dāng)石油泄漏時(shí),泄漏口周圍的空氣可能比石油更冷,因此,泄漏會(huì)導(dǎo)致石油的溫度下降。
為了進(jìn)一步探究泄漏對(duì)整個(gè)管道的沿程流動(dòng)參數(shù)的影響,使用Pipeline Studio軟件對(duì)石油泄漏管道進(jìn)行仿真,仿真數(shù)據(jù)與之前一樣,來(lái)自國(guó)內(nèi)某段油氣管道。泄漏前后管道內(nèi)油液壓力、速度、溫度和密度變化如圖5所示。
圖5 泄漏前后管道沿程流動(dòng)參數(shù)變化曲線
通過(guò)觀察圖5泄漏前后壓力、溫度、密度和速度的變化,不難發(fā)現(xiàn),在輸入端壓力、溫度固定,輸出端流量固定的條件下,當(dāng)管道未發(fā)生泄漏時(shí),油品經(jīng)過(guò)管道流動(dòng)后,壓力、溫度和速度逐漸下降,但由于設(shè)置了保溫層,溫度下降幅度不大。隨著溫度的降低,油品的密度逐漸增大。
當(dāng)管道發(fā)生泄漏后,管道沿程的壓力、溫度和密度的變化趨勢(shì)仍然符合流動(dòng)基本規(guī)律,但是速度會(huì)在泄漏點(diǎn)發(fā)生突降,這是因?yàn)樾孤┦ヒ徊糠至髁?。壓力在泄漏的一瞬間,泄漏點(diǎn)會(huì)發(fā)生突降,然后向管道兩側(cè)傳遞,最后導(dǎo)致整根管道壓力下降幅度增大,末端壓力數(shù)值變小。因?yàn)樾孤┛纵^小,換熱的影響不大,導(dǎo)致溫度和密度基本不受泄漏影響。雖然泄漏前后溫度和密度的變化并不大,但是壓力和速度卻有較大變化,因此,輸入端速度的變化和輸出端壓力的變化將是定位泄漏點(diǎn)的依據(jù)。
上述的泄漏定位原理就是將求解泄漏位置轉(zhuǎn)化為求解最優(yōu)問題,其目標(biāo)函數(shù)是輸入端流速和輸出端壓力的計(jì)算值與實(shí)測(cè)值偏差最小,優(yōu)化參數(shù)就是假設(shè)的泄漏位置和泄漏孔大小,約束條件就是在整個(gè)管道長(zhǎng)度內(nèi)尋優(yōu),并且泄漏孔直徑在0.02~0.06 m,因此,可以得到如下所示的石油管道泄漏定位模型。
(17)
式中: 參數(shù)上帶尖號(hào)的為計(jì)算值,其余為實(shí)測(cè)值。
根據(jù)上一小節(jié)的泄漏定位原理,需要選擇合適的優(yōu)化算法解出管道泄漏孔位置和大小。選擇鯨魚優(yōu)化算法(WOA)求解,并將求解結(jié)果與前人使用過(guò)的遺傳算法(GA)[18]和粒子群算法(PSO)[19]進(jìn)行對(duì)比。
鯨魚優(yōu)化算法(whale optimization algorithm,WOA)是一種啟發(fā)式優(yōu)化算法[20-21],由明Seyedali 于2016年提出。其靈感來(lái)源于鯨魚群體的行為方式,利用了鯨魚的游動(dòng)和社交行為,尋找全局最優(yōu)解。WOA算法通過(guò)迭代優(yōu)化尋找最優(yōu)解。在初始狀態(tài)下,通過(guò)生成一組隨機(jī)的解初始化種群。每次迭代中,WOA算法通過(guò)調(diào)整種群中的個(gè)體逐步優(yōu)化。優(yōu)化過(guò)程中,算法利用了3種行為模式:搜索、盤旋和攻擊。
搜索行為通過(guò)隨機(jī)行動(dòng)發(fā)現(xiàn)新的解。盤旋行為是一種自適應(yīng)的行為,利用了鯨魚的游動(dòng)方式,使種群能夠在解空間中更好地探索。攻擊行為是一種探索局部最優(yōu)解的方式,通過(guò)與種群中最優(yōu)個(gè)體進(jìn)行比較更新個(gè)體的位置。WOA算法通過(guò)自適應(yīng)機(jī)制和隨機(jī)搜索,能夠在優(yōu)化過(guò)程中快速收斂到最優(yōu)解,而且能夠克服陷入局部最優(yōu)解的不足。
以Pipeline Studio仿真結(jié)果[22]作為驗(yàn)證不同優(yōu)化算法差異的依據(jù),Pipeline Studio管道的輸入端和輸出端的壓力和流量變化如圖6所示,該圖描述了600 s內(nèi)起點(diǎn)和終點(diǎn)的壓力和流量的變化趨勢(shì),其中在第50 s該管道發(fā)生泄漏,泄漏點(diǎn)在距離輸入端65 km處,泄漏一直持續(xù)到600 s,且泄漏孔直徑保持為0.02 m不發(fā)生變化。
圖6 輸入端和輸出端的壓力和流量變化曲線
從圖6可以看出,在第50 s發(fā)生泄漏后輸出壓力輸入流量立刻發(fā)生相應(yīng)變化,即輸出壓力下降,輸入流量上升,而且由于泄漏點(diǎn)離輸出端更近,輸出壓力的變化要早于輸入流量;此外,輸入壓力和輸出流量都沒有任何變化,這是由于邊界條件設(shè)定為固定值。當(dāng)時(shí)間到達(dá)450 s時(shí),泄漏工況基本穩(wěn)定,由于沒有實(shí)際生產(chǎn)生活中的噪聲影響,所以各種曲線都比較平滑。
當(dāng)泄漏工況達(dá)到穩(wěn)定狀態(tài)后,需記錄下每秒的數(shù)據(jù),并將每10 s的數(shù)據(jù)歸為一組,連續(xù)記錄600 s,總共將獲得15組數(shù)據(jù)。隨后,將每組數(shù)據(jù)代入到石油管道泄漏定位模型中,采用多種不同的優(yōu)化算法進(jìn)行求解,以得到各種優(yōu)化算法的計(jì)算速度,其結(jié)果如圖7所示。
圖7 不同算法定位時(shí)間對(duì)比曲線
根據(jù)圖7所示,這3種算法都能實(shí)現(xiàn)泄漏管道定位的目的,但是計(jì)算時(shí)間較長(zhǎng),這是因?yàn)槊恳淮味夹枰獙⒁环N泄漏工況迭代計(jì)算到穩(wěn)定狀態(tài)才開始對(duì)目標(biāo)函數(shù)尋優(yōu)。就計(jì)算時(shí)間而言,3種算法中粒子群算法定位時(shí)間最長(zhǎng);鯨魚算法定位時(shí)間最短;遺傳算法雖然定位時(shí)間與鯨魚算法相差不大,但是容易陷入局部最優(yōu)解當(dāng)中,增加計(jì)算時(shí)間。
在之前算法定位時(shí)間的對(duì)比上,雖然3種算法的表現(xiàn)各不相同,但都能實(shí)現(xiàn)對(duì)泄漏管道定位的目的,為避免偶然性,本節(jié)將在不同輸入壓力和輸出流量情況下,驗(yàn)證3種算法的可行性,并探究壓力和流量變化對(duì)算法定位的影響。
通過(guò)改變輸入壓力和輸出流量,模擬不同的泄漏工況,分析不同工況下各種優(yōu)化算法的定位穩(wěn)定性。將輸入壓力設(shè)置為6~7.5 MPa,輸出流量設(shè)置為3 550~3 650 m3/h,泄漏點(diǎn)位置在距輸入端65 km處,泄漏孔直徑為0.02 m,采取和研究定位時(shí)間一樣的方式獲取數(shù)據(jù),得到不同泄漏工況下的定位結(jié)果,見表1。
表1 不同輸入壓力和輸出流量下3種算法的泄漏定位結(jié)果
根據(jù)表1的內(nèi)容,將不同優(yōu)化算法在不同輸入壓力和輸出流量工況下的定位誤差和泄漏孔徑預(yù)測(cè)誤差作圖,如圖8和圖9。
圖8 3種算法在定位上的誤差曲線
圖9 3種算法在孔徑上的誤差曲線
不同的壓力和流量工況下,3種算法都能實(shí)現(xiàn)對(duì)泄漏點(diǎn)的定位和對(duì)泄漏孔徑的預(yù)測(cè)。但是隨著壓力和流量的增大,3種算法并沒明顯的變化趨勢(shì),表明壓力和流量的變化對(duì)算法的影響不大。將3種算法對(duì)比發(fā)現(xiàn):就3種算法對(duì)泄漏孔的定位來(lái)說(shuō),鯨魚優(yōu)化算法的定位誤差平均為0.37%,遺傳算法定位誤差平均為0.43%,粒子群算法定位誤差平均為0.49%,雖然平均誤差看上去相差不大,但是鯨魚優(yōu)化算法更加穩(wěn)定,適用性更好;從3種算法對(duì)泄漏孔大小預(yù)測(cè)來(lái)說(shuō),3種算法的誤差波動(dòng)都比較大,但是鯨魚優(yōu)化算法表現(xiàn)得相對(duì)較好,誤差更小一些。
只是在改變輸入壓力和輸出流量情況下驗(yàn)證3種算法的可行性明顯是片面的。在實(shí)際情況下,管道泄漏的位置和孔徑存在很大的不確定性,為了更符合實(shí)際,本節(jié)將在不同泄漏位置和孔徑情況下,驗(yàn)證3種算法的可行性,并探究泄漏位置和孔徑大小對(duì)算法定位的影響。
通過(guò)改變泄漏位置和泄漏孔徑大小,形成不同的泄漏工況,分析不同工況下各種優(yōu)化算法的定位精度。將泄漏位置分別設(shè)置為20、35、50、65 km,泄漏孔直徑設(shè)置為20、30、40、50、60 mm,采取和研究定位時(shí)間一樣的方式獲取數(shù)據(jù),得到不同泄漏工況下的定位結(jié)果,見表2。
表2 不同泄漏位置和泄漏孔徑下3種算法的泄漏定位結(jié)果
根據(jù)表2的內(nèi)容,將不同優(yōu)化算法在不同輸入壓力和輸出流量工況下的定位誤差和泄漏孔徑預(yù)測(cè)誤差繪制成圖10和圖11。
圖10 3種算法在定位上的誤差曲線
如圖10和圖11所示,不同的泄漏孔位置和泄漏孔徑下,3種算法都能實(shí)現(xiàn)對(duì)泄漏點(diǎn)的定位和對(duì)泄漏孔徑的預(yù)測(cè),并且隨著泄漏位置和泄漏孔徑的增加,可明顯發(fā)現(xiàn)3種算法對(duì)泄漏位置的定位更加準(zhǔn)確。
雖然3種算法都實(shí)現(xiàn)了定位的目標(biāo),但是精度各不相同。就3種算法對(duì)泄漏孔的定位來(lái)說(shuō),粒子群優(yōu)化算法的誤差相對(duì)小一些,但是與鯨魚優(yōu)化算法差距并不大,相比之下,遺傳算法就顯得不夠優(yōu)秀;從3種算法對(duì)泄漏孔大小預(yù)測(cè)來(lái)說(shuō),3種算法的誤差波動(dòng)都比較大,但是鯨魚優(yōu)化算法平均誤差在7%左右,而遺傳算法和粒子群算法的平均誤差都在10%以上,綜合考慮可知,鯨魚優(yōu)化算法的精確性更好。
1) 切應(yīng)力修正過(guò)后的石油管道模型是準(zhǔn)確的。在同樣泄漏工況下,修正后管道模型與pipeline studio軟件計(jì)算的沿程壓力最大相對(duì)誤差為1.02%,沿程流速最大相對(duì)誤差為3.84%。
2) WOA、GA、PSO算法都能求解石油管道泄漏模型,但是WOA算法計(jì)算的時(shí)間更快。3種算法隨著管道輸入端壓力和輸出端流量的變化沒有明顯的改變,但是隨著泄漏位置的改變和泄漏孔徑的增大,定位更加準(zhǔn)確。
3) 隨著管道輸入端壓力和輸出端流量的變化,WOA算法的定位誤差平均在0.37%左右,GA算法的定位誤差平均在0.43%左右,PSO算法的定位誤差平均在0.49%左右;并且WOA對(duì)泄漏孔徑的預(yù)測(cè)更穩(wěn)定。
4) 隨著泄漏位置的改變和泄漏孔徑的增大,PSO算法的定位誤差相對(duì)更小,但是WOA算法與它的差距并不大;并且對(duì)于孔徑預(yù)測(cè),WOA算法平均誤差在7%左右,而GA算法和PSO算法的平均誤差都在10%以上。