董峰 李昭 李凌涵 張淑美
氣液兩相流是一種氣體和液體同時存在且具有明確分界面的流體流動形式,廣泛存在于自然界及能源、動力、石油、化工、冶金、制藥等眾多工業領域中,具有流動結構復雜、過程狀態波動、過程變量多及狀態的變化具有隨機性、非線性等特點.由于氣液兩相流流動狀態隨環境、各分相流體性質及相間的相互作用而改變,其過程參數的檢測和流動狀態的監測相比單相流具有更大的挑戰[1-2].如能及時、準確地獲得流動過程狀態信息,對流動過程的產生、發展及轉化進行分析,實現流動狀態的表征和描述,可為進一步實現氣液兩相流流動過程的控制提供關鍵參數,對促進多相流問題研究的進一步發展,保障實際生產過程的安全穩定運行具有重要意義.
在氣液兩相流流動過程中,流速、相含率、壓力或壓力降等過程參數及相介質的分布反映了其流動過程特性.其中,相含率可以通過電學法[3-4]、射線法[5]、微波法[6]、快關閥法[7]等獲取;流速可以通過差壓法[8]、超聲多普勒方法[9]等測量;截面陣列式電阻傳感器可以獲取流體不同位置處的相分布信息.在獲取相關流動參數及相分布特性等基礎上,采用小波變換[10]、經驗模態分解、希爾伯特-黃變換[11]等時頻分析方法,可以進一步揭示氣液兩相流流動過程狀態的演化特性,為流動狀態的判別和過程監控提供了有益的借鑒.
通過多傳感器獲取同類傳感器的多模式信息以及異類傳感器的多檢測角度信息,可以更加完整地表征氣液兩相流的流動狀態[1].根據不同傳感器特點進行原始信號預處理,獲得的多維時間序列蘊含了流動狀態的變化特性.通過歷史數據建立狀態監測模型,獲取實時數據進行狀態監測,為實現氣液兩相流復雜流動狀態的監測和分析提供了可行性方法.
基于數據驅動的多元統計監測方法如主成分分析(Principal component analysis,PCA)[12]、獨立成分分析(Independent component analysis,ICA)[13]、偏最小二乘(Partial least squares,PLS)[14]等可以簡化數據結構,解決多變量的耦合問題,在實際過程得到廣泛應用[15-16].在氣液兩相流研究領域,PCA、ICA 等方法多用于流型識別和流體測量中傳感器信號的特征提取.Shaban等[17]利用歸一化后的差壓信號建立PCA-ICA 模型,實現了垂直氣液兩相流流量測量;李凱鋒等[18]利用PCA 提取電導率信息主成分,結合K-均值聚類算法實現了氣水兩相流流型辨識;Dong等[19]對電導環信號的多域特征采用等距特征映射、PCA 等方法得到三維向量實現了流型樣本間結構的可視化.
PCA 作為實際工業過程最為成熟和常用的過程監測方法,可以從相互關聯的高維數據中提取有用信息,實現數據的有效降維.Li等[20]將PCA 方法和支持向量數據描述結合起來用于制冷系統的故障診斷;Aouabdi等[21]結合多尺度熵和PCA,通過分析電機電流信號監測齒輪退化.PCA 假設過程變量是線性的且在時間上靜態獨立,僅提取了變量間的靜態互相關,但實際的工業過程數據往往具有非線性、動態性等特點[22-23].Ku等[24]通過增加遲滯數據對數據矩陣進行動態拓展,提出了一種動態主成分分析方法(Dynamic PCA,DPCA).DPCA 提取變量在時間上的自相關和動態互相關關系,但仍未考慮到過程數據的非線性特性.基于以上問題,Choi等[25]提出了基于動態核主成分分析(Dynamic kernel PCA,DKPCA)的非線性動態過程監測方法,Zhang等[26]將特征向量選擇與DKPCA 相結合提高了故障檢測效率.
DKPCA 在工業領域中主要用于監測正常工況是否出現異常故障,時滯數據擴展適合描述過程的動態特性,核函數有利于捕獲過程的非線性特性,其方法本身同時考慮了過程數據的動態性和非線性,符合氣液兩相流流動過程特點.在應用于氣液兩相流流動狀態監測時,更應注重對流動狀態的描述和分析[27-28].
為實現氣液兩相流動過程狀態的監測,以實際流動過程中出現頻率較高的泡狀流、塞狀流和彈狀流3 種典型流動狀態為研究對象,利用不同流動狀態下多傳感器測試數據所呈現的流動特性差異性,提出一種針對氣液兩相流狀態監測的多模態動態核主成分分析方法(Multiple dynamic kernel PCA,MDKPCA).通過對多種流動狀態采用歷史數據建立多模態監測模型,實現對典型流動狀態的判別及過渡狀態的監測.
DKPCA 將過程的動態特性和非線性特性融入模型,將輸入空間X∈Rn×m拓展到具有l個時滯的增廣矩陣X(l)以提取過程數據的動態特性[29],即

構造從輸入空間到特征空間的非線性映射,以提取過程數據的非線性特性[30],在線性特征空間的協方差可以表示為

式中,Φ(·)表示將輸入空間的輸入向量映射到特征空間F的非線性映射函數.
通過求取特征空間的特征值可以對應得到特征空間的線性主元

式中,λ≥0,ν∈F{0},將λ值降序排列,最大的λ對應的ν為特征空間F中的第一個主元,最小的λ對應的ν為最后一個主元.
輸入向量x(l)的主元成分t通過在特征空間F將Φ(·)映射到νk(k=1,2,···,p)而被提取

為避免直接計算非線性映射,引入點積形式的核函數k(x,y)=〈Φ(x),Φ(y)〉.
氣液兩相流中不同的典型流動狀態具有不同的流動特性,體現在過程參數的差異.因此,可采取多模態建模[31],即分別對每一種流動狀態利用歷史數據采用DKPCA 建模,根據實時數據進行狀態監測.以實際流動過程中出現頻率較高的泡狀流、塞狀流和彈狀流3 種典型流動狀態為研究對象,MDKPCA 建模步驟為:

4)將步驟2)中的徑向基核函數k(x,y)代入式(4)中,分別提取3 種模態數據的非線性主成分:

5)計算T2和SPE 統計量及統計量控制限,統計量計算式為

式中,t為保留p個主元的主元得分向量,對角矩陣S=diag{λ1,···,λp}由建模數據協方差矩陣與保留主元數相對應的p個特征值所構成,表示保留所有主元的得分向量的平方和,表示保留p個主元的得分向量的平方和.
根據方差累計貢獻率(Cumulative percent variance,CPV)選擇主元數量.經檢驗,氣液兩相流的過程數據不符合高斯分布,因此選用核密度估計法(Kernel density estimation,KDE)計算T2及SPE 統計量控制限和SPEq,lim.二者的置信度水平均選擇為α=0.99.
氣液兩相流流動過程作為動態系統,狀態與狀態間并非完全獨立,當前的狀態很大程度上依賴于過去的狀態.單一時刻的過程參數在一定程度上并不能準確表征流動狀態的特性:對于塞狀流流動狀態,一個完整的流動周期應包括由細小氣泡聚集而成的長塞型氣泡以及緊隨其后的細小氣泡,二者交替出現,間斷地沿管壁頂部流動;對于彈狀流流動狀態,在一個流動周期內,含大量液體和細小氣泡的液彈和含有大量氣體的氣彈交替出現.當采樣點取到塞狀流的細小氣泡段及彈狀流的液彈段,并不能表征其完整的流動過程特性.因此,在監測時應考慮到之前時刻的數據值,采用滑動窗(Moving window)技術,將當前采樣時刻數據與前l個時刻的采樣數據構成一個窗口寬度為l+1 的新數據樣本.
監測步驟為:
1)獲取新的實時數據xnew,按照所建立的不同狀態模型的均值和方差進行標準化處理,得到經標準化處理后的xq,new.
2)將xq,new與前l個標準化處理后的數據納入滑動窗,作為當前時刻采集到的數據樣本new=[xq,l,xq,l-1,···,xq,new]T.
3)利用建模步驟3)中的Kq和1N,計算核向量kt∈R1×N并進行均值中心化處理:

4)利用建模步驟4)的式(6),令xq(l)=xq,new,提取非線性主成分tq,new.
如圖1 所示為氣液兩相流流動狀態MDKPCA建模與監測原理圖.

圖1 氣液兩相流MDKPCA 建模與監測原理圖Fig.1 MDKPCA modeling and monitoring schematic diagram of gas-liquid two-phase flow
MDKPCA 方法用于氣液兩相流狀態監測的可行性和有效性通過實驗進行驗證.實驗產生不同工況下的多種流動狀態,通過多傳感器采集流動過程數據并進行預處理,進而進行MDKPCA 建模.
氣液兩相流實驗于氣液兩相流水平環管實驗裝置上完成,實驗裝置如圖2 所示.

圖2 氣液兩相流水平環管實驗裝置Fig.2 Experimental apparatus of horizontal loop for gas-liquid two-phase flow
實驗所用的管道為內徑50 mm 的不銹鋼管,流體分別為自來水(密度為998 kg/m3,動力粘度為1.01×10-3Pa · s)和干燥空氣(密度為1.2 kg/m3),在管道入口處安裝的混合器進行混合.流體進入混合器之前,使用標準單相流量計計量各相流量,由入口至出口距離為16.6 m.為獲得充分發展的流動狀態,多傳感器測試管段安裝在位于入口處下游約12 m.測試管段為透明有機玻璃管,可觀察、并采用高速攝像機記錄流動過程狀態.流動過程的工況條件由測試管段上的壓力和溫度儀表進行監測和記錄.
在氣液兩相流流動過程中,管道內各相流速、相分布、相含率及壓力、溫度等過程參數可表現流動狀態的特性.為了多方面獲取流動過程信息,選擇用于數據獲取的多傳感器有:截面陣列式電阻、連續波超聲多普勒、電容、電導環以及壓力計.測試管段多傳感器布置結構如圖3 所示.各傳感器所獲數據表征的流體特性信息及相關的過程參數如表1所示.

圖3 測試管段多傳感器結構Fig.3 Structure of multiple sensors in test section

表1 氣液兩相流測量多傳感器Table 1 Multiple sensors for gas-liquid two-phase flow measurement
實驗包含7 組、共42 個實驗條件測試點.每組水流量固定不變,通過調節空氣流量由低至高變化,產生不同的流動狀態,包括泡狀流、塞狀流、彈狀流3 種典型流動狀態和從泡狀流到彈狀流的過渡狀態.水的流量依次設置為2,5,7,9,11,13和14 m3/h,標況(0 攝氏度,1 個大氣壓)下空氣流量為1.5,2.5,5,7,19和36 m3/h.圖4 所示為實驗點在Mandhane 氣液兩相流流型圖中的分布[32].圖中實線為流型轉換邊界,不同圖標點表示實驗中出現的典型流動狀態.

圖4 氣液兩相流實驗點分布Fig.4 Distribution of experimental points for gas-liquid two-phase flow
1)截面陣列式電阻
截面陣列式電阻傳感器的工作原理為通過激勵電極在測量區域建立敏感場,利用流體流動狀態改變引起敏感場域內的電導率改變,導致測量邊界電壓發生變化.16 個電極沿著管壁等距離排列,采用相鄰激勵、相鄰測量模式,采樣頻率為140 幀/s,每幀包含208 個邊界電壓測量值,對測量邊界電壓數據進行預處理[33]

式中,Vij0表示管道內充滿水時第i個激勵電極下的第j個邊界電壓測量值,Vij表示兩相流動時第i個激勵電極下的第j個邊界電壓測量值,VRi表示兩相流動時第i個激勵電極下13 個邊界電壓的平均值.
由于氣液兩相流流動結構復雜,截面陣列式電阻的測量數據反映了流體流動過程的二維介質分布信息,因此1 幀的數據并不能完整地表征流動過程的狀態.根據測試系統中其他傳感器的測試頻率,需進行測試數據的時間配準.采用每次移動1 幀的滑動窗方法,對連續42幀VRi數據的幅值進行特征提取,并實現數據降維.提取的特征值包括均值、方差SD、偏度SK和峭度KI等用于表征流動氣、液相介質的分布信息

式中,xi為時間序列VRi中的采樣點,為xi的樣本均值,n為時間序列VRi的采樣點數量.
2)連續波超聲多普勒
連續波超聲多普勒傳感器由一對壓電陶瓷超聲波換能器構成.根據管徑及流速范圍,選用直徑9 mm 的壓電晶片,安裝傾斜角為55°.超聲波發射器以1 MHz 的頻率向被測流體中發射連續聲波,超聲波接收器接收被測流體中經離散相液滴調制的聲波.根據超聲多普勒效應,多普勒頻移大小直接取決于離散相液滴的流速.連續波超聲多普勒采樣頻率為50 kHz,采集的電壓數據為時間序列,采用短時傅里葉變換(Short-time Fourier transform,STFT)將時域信號轉變為頻域信號,并求取多普勒速度[34]

式中,Sd(fd)為多普勒頻移的能量譜,fd為多普勒頻移組分,f0為超聲換能器發出的超聲波頻率,c為超聲波在流體中的縱波波速,θ為超聲波聲束方向與水平方向的夾角.
3)電導環
電導環傳感器由軸向排列的6 個環形金屬電極組成,采樣頻率為2 kHz,其中,E、G 為激勵電極,M 為測量電極,相鄰電極之間的距離為:E -M1和M4 -G 為50 mm,M1 -M2和M3 -M4 為20 mm,M2 -M3 為60 mm.采用電流激勵電壓測量的方式,通過激勵電極對E-G,向管道內施加20 kHz的方波激勵電流形成電場.連續相導電時(水為連續相),可以通過測量電極對之間的電壓表征連續相(水)含率信息[34]

式中,Vw為滿水測量值,Vmeas為實際測量值,Vn為歸一化電壓值.
4)電容
電容傳感器由一對凹面金屬極板C1和C2 構成,極板間混合流體的介電常數由各相的介電常數和相分布決定,因此測量電容與相含率存在對應關系.當連續相不導電時(氣為連續相),可通過電容值計算離散相(水)含率.為消除管壁對檢測結果帶來的影響,采用相對電容變化量RCD表征離散相(水)含率信息[35]

式中,Cm為流體測量電容值,Ct為管道充滿導電相時的電容值,CI為管道充滿非導電相時的電容值,Vm為流體測量電壓值,Vt為管道充滿導電相時的電壓值,VI為管道充滿非導電相時的電壓值.
壓力計測量管道內壓力,主要反映了流動過程的工況信息.由于環境溫度變化幅度很小且對氣液兩相流流動狀態影響較小,故未考慮溫度參數.
各傳感器每個原始數據集的采樣時間均為10 s,但由于采樣頻率不同,因此通過重采樣對齊多傳感器樣本數據.經重采樣后,每組的數據為Y(1 000×8)矩陣.以此為基礎,對泡狀流、塞狀流、彈狀流3 種典型流動狀態進行MDKPCA 建模,并分別獲取各典型流動狀態下的T2 及SPE 統計量控制限.考慮到實際過程以及實驗所獲取的各流動狀態的采樣數量,采用的建模數據:泡狀流為5 000 個采樣點,塞狀流為4 000 個采樣點,彈狀流為10 000 個采樣點.由于氣液兩相流流動過程的狀態需要一段時間內結構的變化進行描述,因此在分析過程狀態的動態特性時,需采用一定時長的測試數據.根據實際氣液兩相流流動過程狀態變化的特點,取采樣的時滯l=29.即以時間序列長度為30 (對應時間長度為0.3 s)的采樣點為一個樣本點用于建模和測試.
以氣液兩相流在實際流動過程中最為常見的泡狀流、塞狀流、彈狀流等典型流動狀態為研究對象,3 種典型流動狀態如圖5 所示.泡狀流中,氣泡非均勻分散于管道的中上部;塞狀流中,部分細小氣泡聚集成長氣塞,氣塞尾部伴隨細小氣泡在管道頂部流動;彈狀流中,夾帶氣泡的液彈與較大速度的沖擊性氣彈快速交替.

圖5 3 種典型流動狀態Fig.5 Three typical flow states
分別采用多模態的主成分分析(Multiple PCA,MPCA)、多模態的動態主成分分析方法(Multiple dynamic PCA,MDPCA)和MDKPCA 方法進行狀態判別,結果分別如圖6~8 所示.
圖6 中,由于采用的MPCA 方法僅提取變量間的靜態互相關信息,組成數據樣本的時滯l=0,樣本數據為:1~1 000 采樣點為泡狀流狀態,1 001~2 000 采樣點為塞狀流狀態,2 001~3 000 采樣點為彈狀流狀態.MPCA 狀態判別結果表明:泡狀流在3 種監測模型下除個別點外均未超限,即泡狀流符合泡狀流、塞狀流及彈狀流3 種流動狀態模型的特性;此時,會將泡狀流誤判為塞狀流或者彈狀流.這一現象符合3 種流動狀態的流動過程規律:塞狀流和彈狀流的完整周期包括細小氣泡段.MPCA 方法只針對單一時刻的狀態進行獨立的線性分析,無法描述流體流動過程中前后狀態間的關聯性和狀態變化的非線性,導致泡狀流的誤判.而塞狀流和彈狀流整體上判別效果好于泡狀流,可以比較直觀地區分,但在非塞狀流和彈狀流對應的監測模型下,仍有大量采樣點在控制限下,存在較高的誤判率.

圖6 MPCA 監測模型典型狀態判別Fig.6 MPCA-based state identification for typical states
圖7 中,采用MDPCA 方法進行3 種流動狀態的判別時,組成數據樣本的時滯l=29,樣本數據為:1~971 樣本點為泡狀流狀態,972~1 942 樣本點為塞狀流狀態,1 943~2 913 樣本點為彈狀流狀態.MDPCA 狀態判別結果表明:泡狀流在3 種監測模型下均未超限,但相比于MPCA 方法,泡狀流樣本點在非泡狀流監測模型下,均比模型對應的塞狀流或彈狀流數據整體更接近控制限.

圖7 MDPCA 監測模型典型狀態判別Fig.7 MDPCA-based state identification for typical states
圖8 中,采用MDKPCA 方法對3 種流動狀態進行判別時,組成數據樣本的時滯l=29,樣本數據為:1~971 樣本點為泡狀流狀態,972~1 942 樣本點為塞狀流狀態,1 943~2 913 樣本點為彈狀流狀態.MDKPCA 狀態判別結果表明:泡狀流在塞狀流監測模型下,SPE 統計量完全超限;在彈狀流監測模型下,雖仍有部分樣本點低于控制限下,但相比于彈狀流數據,T2和SPE 統計量整體更接近控制限.塞狀流和彈狀流數據的誤判率進一步降低.相比于前兩種方法,3 種典型流動狀態的判別效果有了很大的提升.

圖8 MDKPCA 監測模型典型狀態判別Fig.8 MDKPCA-based state identification for typical states
為進一步討論不同典型流動狀態對MPCA、MDPCA和MDKPCA 等3 種方法所建立監測模型判別效果的影響,對T2和SPE 統計量進一步分析.定義模型對各典型流動狀態的符合度為:觀測的典型流動狀態被判別為所用模型對應流動狀態的比率.定義判別正確率β為:3 種監測模型對觀測的典型流動狀態均正確判別的比率.
分別采用MPCA、MDPCA和MDKPCA 方法得到各監測模型對所觀測到的3 種典型流動狀態的判別結果如表2 所示.
表2 所示數據中,MDPCA 方法對塞狀流和彈狀流的判別效果明顯優于MPCA 方法,將塞狀流和彈狀流在非對應監測模型下誤判為其他流動狀態的比例明顯降低,判別正確率明顯提高.表明了MDPCA 方法通過提取狀態變化的動態自相關和互相關信息,把流體流動過程近似為一種線性的動態過程處理后,對過程數據中包含的動態特性敏感性增強.但MDPCA 方法難以將泡狀流和其他兩種流動狀態的流動特性區分,觀測的泡狀流狀態在塞狀流和彈狀流監測模型下仍存在100%和76.93%的誤判率.分析原因,主要是MDPCA 方法未考慮流動過程的非線性變化及測試數據中耦合的這種非線性變化特點.

表2 觀測流動狀態下3 種方法所建立監測模型判別結果(%)Table 2 Identification results of monitoring models in three methods under observation of flow states (%)
MDKPCA 方法在MDPCA 方法的基礎上,進一步將測試數據中耦合的非線性變化特點,通過核方法映射到高維空間進行線性化處理,提取塞狀流和彈狀流在流動過程中液滴和氣泡不斷聚并、破碎的動力學演化過程的細節特征,用于區分不具有上述演化過程特征的泡狀流狀態.因此,MDKPCA方法與MPCA和MDPCA 方法相比,建立的監測模型對模型非對應的典型流動狀態誤判率明顯降低,3 種典型流動狀態判別正確率明顯提高.MDKPCA 方法可以更好地描述氣液兩相流不同流動狀態的流動過程特性,進行狀態監測時效果更佳.
在氣液兩相流流動過程中,從一種穩定的流動狀態演化到另一種穩定的流動狀態之間存在過渡狀態,以穩定的泡狀流過渡到穩定的彈狀流這一完整的流動過程為研究對象,由高速攝像機記錄觀測到的過程狀態如圖9 所示.分別采用MPCA、MDPCA、MDKPCA 方法進行狀態監測,監測結果分別如圖10~12 所示.

圖9 泡狀流過渡到彈狀流的過程狀態Fig.9 Transition from bubble flow to slug flow
圖10 中,采用MPCA 方法進行狀態監測,1~2 000 個采樣點為穩定的泡狀流,2 001~4 000 個采樣點為過渡狀態,4 001~6 000 個采樣點為穩定的彈狀流.MPCA 監測結果表明:在泡狀流監測模型下,只有T2統計量可以體現從穩定的泡狀流過渡到穩定的彈狀流狀態轉變的過程;在彈狀流監測模型下,兩個統計量值在整個過程中均位于對應的控制限下,完全無法區分流動狀態和實現狀態監測.

圖10 MPCA 監測模型過渡狀態監測Fig.10 MPCA-based monitoring for transitions
圖11 中,采用MDPCA 方法進行狀態監測,1~1 942 個樣本點為穩定的泡狀流,1 943~3 884 個樣本點為過渡狀態,3 885~5 826 個樣本點為穩定的彈狀流.MDPCA 監測結果表明:在泡狀流監測模型下,兩個統計量均可以體現從穩定的泡狀流過渡到穩定的彈狀流狀態轉變的過程;在彈狀流監測模型下,雖然兩個統計量都可以反映整個過渡過程的變化趨勢,但不能準確區分泡狀流穩定狀態.

圖11 MDPCA 監測模型過渡狀態監測Fig.11 MDPCA-based monitoring for transitions
圖12 中,采用MDKPCA 方法進行狀態監測,1~1 942 個樣本點為穩定的泡狀流,1 943~3 884 個樣本點為過渡狀態,3 885~5 826 個樣本點為穩定的彈狀流.MDKPCA 監測結果表明:在泡狀流監測模型下,兩個統計量均直觀地體現了整個流動過程的狀態變化;在彈狀流監測模型下,SPE 統計量可以準確區分泡狀流,實現對穩定和過渡狀態的準確監測.當處于過渡狀態時,統計量值的移動和波動反映了在泡狀流過渡到彈狀流過程中,氣相流量不斷增大,氣相和液相呈現出劇烈、復雜的相互作用,揭示了氣液流動過程的轉化并非是瞬態的,而是過去狀態的不斷累積和演化.

圖12 MDKPCA 監測模型過渡狀態監測Fig.12 MDKPCA-based monitoring for transitions
針對水平氣液兩相流流動狀態監測問題,從其狀態的變化具有時變性、非線性、隨機性等特點出發,在多傳感器獲取流動狀態測試數據的基礎上,通過采用動態自相關及互相關方法提取過程的動態特性,核方法提取非線性特性,結合主成分分析的MDKPCA 方法,建立不同典型流動狀態的監測模型,并以T2和SPE 統計量作為監測指標,采用滑動窗技術保障氣液兩相流流動過程信息的連續性,結合所建立的多模態模型對不同的典型流動狀態進行判別,并進一步實現流動過渡狀態監測.
與采用MPCA和MDPCA 方法相比,MDKPCA方法在前者的基礎上,通過進一步提取流動過程的動態和非線性特征,能夠更全面地描述不同流動狀態的過程特性.氣液兩相流實驗裝置測試數據的處理結果驗證了MDKPCA 方法對典型流動狀態判別的準確性和對過渡狀態監測的有效性.