小通道氮氣-水兩相流三譜切片波動特性及其流型表征*

李洪偉 周云龍 王世勇 孫斌

1)(東北電力大學能源與動力工程學院,吉林 132012)

2)(撫順礦業中機熱電有限責任公司,撫順 113001)

(2013年1月22日收到;2013年4月7日收到修改稿)

1 引言

小通道兩相流廣泛地存在于石油、化工和能源等現代工業領域.與常規通道相比,微尺度通道的尺寸減小2—3數量級,使其兩相流動特性表現出與常規通道尺度下的流動特性有很大區別,這主要是因為在微小通道內兩相間的表面張力和潤濕性等處于主導地位.由于氣液兩相間內在的復雜性,有關微小通道內氣液兩相流動特性的研究僅僅處于初級階段,但是因其在化學、化工和能源方面得到廣泛的應用,所以一直以來都受到人們的廣泛關注[1,2].

兩相流是一個復雜的非線性動力學系統,自20世紀90年代以來,基于混沌、分形、復雜網絡以及時頻域分析等方法進行流型識別的研究日趨增多[3],文獻[4—7]采用多個混沌參數指標、有限穿越可視圖以及相空間多元圖重心軌跡等方法分析了油水和氣液兩相流型電導信號的動力學特性,在揭示兩相流型流動機理方面取得了一定進展.Sun等[8,9]應用自適應最優核時頻分析方法對水平文丘里管內氣液兩相流型與流動特性進行了識別與分析.Du等[10]同樣應用自適應最優核時頻表征方法對氣液兩相流電導信號進行了分析,通過提取總能量和時頻熵兩個特征值準確區分出不同流型.Manfredo[11]基于氣液兩相流圖像信號相密度時頻分布特征對其流動特性進行了表征.Ommen等[12]對兩相流時頻域分析方法進行了綜述,分析了在表征兩相流流型動力學特性方面的各自優勢與不足.

高階譜估計技術是分析非高斯信號的有力工具,它從更高階概率結構表征隨機信號,可以彌補二階統計量(功率譜)不包含相位信息的缺陷.高階譜有很強的消噪能力,能夠辨識非因果、非最小相位系統,檢測和表征系統的非線性等[13].高階譜中Welch,AR和EV功率譜近年來被廣泛應用,但幾種分析方法或多或少存在一些問題,最為普遍的缺點是抗噪能力不強且喪失了相位信息,這在非線性系統分析中是致命的弱點.而由文獻[14]演化而來的三譜估計可以有效地規避上述問題,其1.5維切片對系統內在動力學特性揭示更加直觀及具體.

目前流型研究大都集中在常規管道和微米級通道,而介于兩種尺度之間的毫米級小通道(尤其是1.0—5.0 mm)的研究成果相對匱乏,相關文獻報道很少.基于此,本文應用基于三譜切片的波動特性理論,針對內徑2.0 mm×0.81 mm矩形垂直玻璃管道內的氣液兩相流流型動力學特性進行了研究.

2 理論基礎

2.1 高階譜

高階譜分析是一種新的信號處理的有效技術,它可以定量描寫信號間的非線性耦合,可以抑制噪聲、保留信號的相位信息,是處理非線性、非高斯信號的有力工具[15,16].

設系統的傳遞函數為H(X),而a(t)為零均值的非高斯白噪聲,由于采樣信號常受到加性高斯噪聲u(t)的干擾,用傳統的時間序列建模無法排除加性的高斯噪聲干擾.設y1(t)為y(t)中去掉確定性信號的零均值的有色隨機信號,利用y1(t)的高階累積量進行建模就可以抑制u(t)的干擾.對(2)式中的高階累積量 cy,k(τ1,τ2,···,τk-1) 進行 k-1 維的離散Fourier變換,得到高階譜的一般定義為

(1)式中,ω為頻率,H(ω)為系統的傳遞函數,H?(ω) 為 H(ω) 的共軛函數,Sy,k(ω1,ω2,···,ωk-1)為高階譜,也被稱為高階累積量譜.當k=2,3,4時,它分別表示功率譜 P(ω)、雙譜B(ω1,ω2)和三譜T(ω1,ω2,ω3).

2.2 自回歸三譜

自回歸(AR)模型,

式中,ai(i=1,2,···,p)為自回歸系數,p為自回歸模型的階數.對于穩定的線性物理過程h(t),考慮到系統為最小相位系統,根據(1)式得到基于AR模型的雙譜和三譜表達式分別為

對y1(t)運用AR模型的參數法估計出模型系數a,結合(4)式將系數a代人(2)式和(3)式,得到AR雙譜和三譜幅值為

2.3 切片譜

二次相位耦合是常見的一種現象.對于這種非線性耦合現象,僅用功率譜很難從根本上解決問題,而高階譜就可以定量地描述二次相位耦合程度.當凍結三譜 TAR(ω1,ω2,ω3)中的一個頻率 ω3,即當ω3=C(const)時,三譜切片就變成了雙譜[17].定義三譜二維切片為

其規一化后的幅值為

用三譜切片譜函數計算二次相位耦合(信號兩種主要振蕩模式相互耦合),存在計算量大、三維圖形不夠直觀等缺點.為了改善圖形顯示不夠直觀的缺點,采用二維的對角切片譜函數表征系統的1二次相位耦合程度.當ω1=ω2時,定義雙譜的12維(對角)切片譜函數為

2.4 波動特征參數

概率密度分布為非正態分布的隨機信號統稱非高斯信號,在工程中通常用偏斜因子S,平坦因子K和波動強度σ/〈p〉三個參數來描述.本文將這三個參數應用到三譜切片的分析中,實現了三譜切片譜的定量分析.

式中x為離散時間序列,N為序列長度,E(·)為括號內均值.波動強度反映了序列幅值變化梯度的強烈程度,偏斜因子是衡量隨機信號的分布偏離對稱分布的歪斜程度,偏斜因子不等于零的信號必定服從非對稱分布;而平坦因子表征統計頻率曲線接近分布中心時的大致狀態.這三個參數都從不同方面反映出了序列的復雜程度.

3 方法評價

為了評價本文提出的方法在時間序列的處理分析中的有效性,本文主要將抗噪性與復雜度表征性作為評價指標,從具有代表性的分形、周期、混沌時間序列著手考察該方法在處理不同類型時間序列所表現出的抗噪與表征能力.

3.1 分形時間序列

本文對Brown分形序列進行分析,考察三譜切片波動方法在處理分形序列時抗噪與辨識方面所表現出的能力.

分形Brown運動(fractional Brownian motion,fbm)是典型的隨機性分形,其方程 f(x)是一個實值隨機方程{：}

其中,x是E維歐空間RE中的一點;F(t)為隨機變量t的分布函數,該隨機變量服從標準正態分布 N(0,σ2);‖Δx‖ 為采樣間隔;H 為 Hurst指數[6].本文采用Matlab中的wfbm函數生成長度分別為10000,5000,1000,500的分形Brown序列各5組,分別取Hurst指數H=0.1,H=0.3,H=0.5,H=0.7,H=0.9,結果示于圖1.

本文對Brown分形序列的分析有兩個目的：一是驗證三譜切片的抗噪性能,二是討論三譜切片計算過程中數據的長度與快速傅里葉變換(FFT)長度對辨識結果的影響.

圖1 分形Brown運動時間序列

圖2 是Brown序列三譜圖,Brown分形序列雖然表面上看雜亂無章,但從圖中還是可以看到Brown分形序列的內在規律性.二次耦合點主要集中在中心部位,接近于高斯分布.圖3是分形Brown序列在不同噪聲強度干擾下的三譜切片圖.從圖3中可以看出,五種強度噪聲干擾下的功率譜線趨勢基本相同,只存在微小的細節波動,再次驗證了三譜切片良好的抗噪性能.

在計算同一種時間序列的三譜切片及波動特征值過程中,我們發現有兩個參數對計算結果的影響比較大,一個是所取的序列長度N,另外一個是在進行FFT變換過程中所選擇的長度nFFT.本文對分形Brown運動時間序列在取不同H值情況下進行了N與nFFT的正交計算,觀察每相鄰兩組H值的三個波動特征參數類間距值,長度為10000時的具體數據如表1所示.

圖2 分形Brown序列三譜圖(H=0.1,N=10000)

圖3 Brown序列不同加噪下三譜圖(H=0.1,N=10000)

從表1中可以很清楚地看到,在沒有噪聲干擾下的分形Brown序列隨著nFFT的減小三個波動特征參數間的差值逐漸減小,當nFFT為64時,每相鄰兩組序列的間距值只有一組不為0,且這一組中也只有波動強度值不為0,平坦因子和偏斜因子也為0.這表明在同一長度N下nFFT值越大,波動特征參數間的差值越大,辨識度越高.對于不同長度的分形Brown序列又有什么規律,在加入噪聲干擾情況下規律是否發生變化,圖4給出了詳細的說明.從圖4中可以看出,在沒有噪聲干擾的情況下不同長度序列在取得最大類間距時的nFFT長度都固定為1024,是一條平行于x軸的直線,而在有噪聲干擾的情況下,呈現出近似反比關系,當N越大nFFT取值越小.

表1 分形Brown序列不同n FFT時三個參數間距值

圖4 Brown序列取得最大類間距值下的N與n FFT曲線

表明在時間序列嚴格符合高斯正態分布的情況下,N與nFFT兩者越大辨識效果越好,在有不規則噪聲干擾從而破壞正態分布結構的情況下,N與nFFT成反比關系.這也為后面流型差壓序列分析時的參數選擇提供了理論依據.

3.2 周期、混沌時間序列

本文為了考察三譜切片及其波動理論適用的廣泛性,對正弦、正弦+噪聲混合、Lorenz以及Lorenz+噪聲混合四種典型信號進行分析計算.幾種時間序列的產生條件如下：

1)正弦信號,y1=sin(x),采樣間隔為π/50;

2)正弦+噪聲,y=y1+p×y2,其中y1為正弦序列,y2為白噪聲序列,p為隨機成分的比例,這里取p=0.2;

3)Lorenz混沌信號,由Lorenz方程

初始條件x=2,y=2,z=20,采用四階龍格-庫塔方法迭代,取變量x為仿真序列[18];

4)Lorenz+噪聲,3)中所得變量x+40 dB噪聲.

圖5 Lorenz與sin信號的加噪前后三譜切片圖

圖5 為抗噪性能檢測,從圖中很容易看出,三譜切片在處理混沌序列時抗噪效果非常好,加噪后與加噪前波動趨勢幾乎如出一轍,而在處理周期序列時受噪聲擾動影響十分明顯,抗噪性能相對較差.而圖6則為四種信號三譜切片的波動特征值,從圖中可以看出,每種信號的三個波動參數值排序為：

Lorenz+噪聲＞Lorenz＞sin+噪聲＞sin,這也正表征出四種信號的內在復雜程度.

圖6 Lorenz與sin信號的加噪前后波動特征值分布

通過對分形序列、周期序列和混沌序列的應用分析,得出結論：三譜切片結合波動特性理論在處理分形和混沌序列時具有較好的抗噪能力,處理周期信號時抗噪能力相對較差.在復雜程度表征方面,無論應用在哪種類型的序列中都具有很好的效果.此外,在計算三譜切片過程中,序列長度N與nFFT的選擇應遵循如下原則：在時間序列嚴格符合高斯正態分布的情況下,N與nFFT兩者越大辨識效果越好,在有不規則噪聲干擾從而破壞正態分布結構的情況下,N與nFFT近似成反比關系.

4 流型辨識及其動力學機理分析

實驗系統如圖7所示.該實驗系統主要包括兩部分,即氣液兩相流動與控制系統、數據采集系統.實驗段采用透光性較好的有機玻璃制成.在通道相距200 mm處取兩個直徑為0.5 mm的取壓小孔,上游小孔距通道入口100 mm,下游小孔距通道出口80 mm,矩形通道截面寬為2.00 mm,縫隙為0.81 mm.

實驗中數據的采集主要有壓差信號、單相氣相流量和單相液相流量.其中：壓差信號采用羅斯蒙特3051C差壓變送器進行測量,測量精度為±0.075%,更新響應時間為100 ms/次,差壓校驗量程為0.5inH2O-2000psi,實驗中采集壓差信號的采樣頻率為256 Hz;液體流量測量采用HW5智能金屬管浮子流量計,此流量計以MCU微處理器為核心,瞬時值和累計值的精度為±0.5%;氣體流量測量采用型號為D-600MD數字式氣體流量計,該流量計的量程范圍為0—5 SLM,精度等級為±1%.實驗中實驗環境壓力為0.1 MPa±0.0005 MPa,實驗環境溫度為25°C±0.5°C,表觀氣速為0.1—30 m/s,表觀液速為0.01—5 m/s.

在圖7中實驗段上采集到四種典型流型的差壓信號,并進行了三譜切片計算,結果如圖8所示.

圖7 氮氣-水兩相流實驗系統

圖8 氮氣-水兩相流典型流型差壓信號、三譜圖、切片圖 (a)毛細泡狀流;(b)分散泡狀流;(c)段塞流;(d)環狀流

從圖8中可以看出,不同流型的三譜圖表現差異較為明顯,這是由于不同流型氣液兩相界面分布不同引起壓降變化不同的結果.相較于三譜3D立體圖的復雜等頻對角切片更加直觀,更能體現不同流型的相分布特點.主要通過對二次耦合譜線的判斷挖掘氣相與液相兩種振蕩頻率與相位的耦合特性.小通道內氮氣-水分布情況最為復雜的是環狀流,其對角切片呈現出頻率峰值分布混亂的特點,主要頻率峰值出現在-π/2,-7π/16,-π/4,π/4以及π/16處,其他很多處也存在微小的耦合現象,說明氣液兩相振蕩并無明顯耦合處,對壓降分布的影響也比較隨機,主要是由于環狀流是由段塞流演變而來,在段塞流中氣流量繼續增加使得氣塞被擊碎,形成不規則振蕩式運動模式,附在小通道表面的液膜在氣相推動下呈現出混沌、分形特征的抖動,而這種抖動的不規則性類似于信號當中加入了白噪聲,其二次相位耦合也相對最為復雜.雖然耦合峰值不穩定,但從圖8(d)可以看到環狀流在不規則的切片譜上還是存在幾處相對較為明顯的峰值.這種不穩定峰值形成的原因有兩種：一方面是實驗系統中氮氣瓶輸出氣體時的波動以及氣體流經三通時分配不均勻形成的波動.而這種波動反映到實驗段中存在一定的延遲,使得實驗段中環狀流的液膜被氣相沖擊產生振蕩存在一定的間歇性,這樣就形成了切片譜中的耦合雙峰值;另一方面是環狀流的波動近似混沌、分形并摻有白噪聲,這種不規則振蕩信號本身存在一定的相位周期性,這種周期性使得環狀流切片譜中出現類似圖3和圖5的幾處峰值.

從圖8(a)和圖8(b)中可以看到毛細泡狀流在-π/4到π/4之間存在四處明顯峰值,而分散泡狀流外側兩個峰值與毛細泡狀流較相近,中間兩個峰值相對較小,這與兩種不同類型的泡狀流氣泡尺寸大小與間隔頻率不同有關,毛細泡狀流中氣泡尺寸分布均勻,間隔頻率具有很好的規律性,所以其兩相振蕩頻率耦合譜線明顯,而分散泡狀流氣泡尺寸分布不如毛細泡狀流均勻.同時兩種泡狀流的運動特征存在明顯的差異.氣液兩相流體以毛細泡狀流型流動時,生成的均勻氣泡與矩形管道的兩壁相貼,氣泡的截面積和小通道的橫截面積近似等同,液相不再連續.隨著氣流量的波動,所生成的氣泡時間間隔大概在平均值的上下波動,這樣氣泡與液塞以幾種周期間歇性運動.由于三譜切片主要考察的是信號本身兩種主要振蕩模式的耦合現象,將高于平均值的時間間隔定位一種周期,那么平均值以下的是另外一種周期,所以從切片圖中可以看到明顯四峰值(兩種周期).而分散泡狀流由于不同尺寸的氣泡均單獨且分散地存在于連續液相之中,不具備毛細泡狀流那樣良好的周期間歇性,但由于氣體流量的波動性,使得微小氣泡在隨著液相的運動過程所出現的量不完全一致,時多時少,正是這種原因使得分散泡狀流的差壓波動信號主振蕩模式被氣泡數量大和小的兩種群體的間歇性出現所主宰,所以其切片譜中的雙峰值現象明顯,由于周期時間較長,所以峰值出現在外側.而靠內側的兩個小的峰值是散落在液相當中小氣泡連續不間斷出現使得差壓信號所產生的微小波動所致,由于周期性不明顯,所以峰值不夠突出.而段塞流由于其特殊的流動結構決定了它的氣液兩相振蕩模式耦合特別明顯(如圖8(c)所示),而且由此可以給出段塞流動的頻率,在π/3左右(一種周期).

圖9給出了四種典型流型差壓信號的三譜切片波動特征值分布情況,通過特征值的分布可以實現對流型的辨識.通過圖9(d)可以發現S與K相結合對段塞流的分類效果要明顯好于其他兩組組合,說明針對周期性較強的時間序列,由于其與高斯序列差距明顯,應用S與K處理分析時會與其他接近分形、混沌的時間序列有較明顯的差距,而在毛細泡狀流、分散泡狀流與環狀流的識別上效果不是很理想,主要原因是由于這三種流型差壓信號波動周期性不強,具有分形與混沌特征,在這種情況下,σ/〈p〉就顯示出了它在處理接近高斯分布序列的優勢,從圖9(b),(c)可以看到,環狀流的分類效果相對還是比較理想,但這三個波動特征參數任意兩兩組合的效果都不是十分完美,始終在毛細泡狀流與分散泡狀流的分類上存在一些瑕疵,總有部分流型樣本混雜在一起,難以準確識別.而將三者結合到一起得到了對四種流型的滿意識別結果(圖9(a)).

5 結論

本文提出了一種基于三譜切片的波動特征分析理論,為了驗證該方法的適用性,將其應用到分形、周期以及混沌時間序列,主要驗證了該方法的抗噪能力與對復雜度表征能力.結果表明該方法在處理分形與混沌時間序列時抗噪能力比較突出,而對于周期序列相對較差,基于三譜切片的波動特征值 (σ/〈p〉,S,K)能夠很好地表征出系統的復雜程度.同時發現在時間序列嚴格符合高斯正態分布的情況下,N與nFFT兩者越大辨識效果越好,而在有不規則噪聲干擾從而破壞正態分布結構的情況下,N與nFFT近似成反比關系.

通過將該方法應用到小通道氮氣-水兩相流差壓信號中,發現四種典型流型三譜切片差異明顯,主要振蕩模式的耦合現象可以揭示出不同流型的內在流動機理,環狀流較泡狀流與段塞流更為復雜,其氣相與液膜相互作用顯示出強烈混沌特征,段塞流由于明顯的周期性使得對角切片分布最為簡單明了.

流型差壓信號的三譜切片結合三個波動特征值實現了對四種典型流型的準確識別.其中S,K對非高斯信號識別效果比較突出,而σ/〈p〉對高斯信號的辨識度更高,將三者結合到一起可以使各自優勢互補,實現小通道內氮氣-水兩相流型的準確識別.如將其應用于其他多相流流型識別領域,以及在其他模式識別領域進行拓展,也將是有益的探索與補充.

[1]Ding GL,Huang D P,Zhang CL 2002 Journal of Refrigeration 2 7(in Chinese)[丁國良,黃東平,張春路2002制冷學報2 7]

[2]Li H W,Zhou Y L,Liu X,Sun B 2012 Acta Phys.Sin.61 030508(in Chinese)[李洪偉,周云龍,劉旭,孫斌2012物理學報61 030508]

[3]Gao Z K,Jin N D 2008 Acta Phys.Sin.57 6909(in Chinese)[高忠科,金寧德2008物理學報57 6909]

[4]Zong Y B,Jin N D,Wang Z Y 2009 Acta Phys.Sin..57 7544(in Chinese)[宗艷波,金寧德,王振亞2009物理學報57 7544]

[5]Dong F,Jin ND,Zong Y B 2008 Acta Phys.Sin.57 6145(in Chinese)[董芳,金寧德,宗艷波2008物理學報57 6145]

[6]Zhou T T,Jin N D,Gao Z K 2012 Acta Phys.Sin.61 030506(in Chinese)[周婷婷,金寧德,高忠科2012物理學報61 030506]

[7]Zhao JY,Jin N D 2012 Acta Phys.Sin.61 094701(in Chinese)[趙俊英,金寧德2012物理學報61 094701]

[8]Sun B,Wang EP,Ding Y 2011 Chinese Journal of Chemical Engineering 19 243

[9]Sun B,Wang E P,Zheng Y J 2011 Acta Phys.Sin.60 014701(in Chinese)[孫斌,王二朋,鄭永軍2011物理學報60 014701]

[10]Du M,Jin N D,Gao Z K 2012 Chem.Engin.Sci.82 144

[11]Manfredo G 2013 Int.J.Multiphase Flow 51 1

[12]van Ommen JR,Sasic S,van der Schaaf J,Gheorghiu S,Johnsson F,Coppens M 2011 Int.J.Multiphase Flow 37 403

[13]Barker RW,Klutke G,Hinich M J1993 ASMETrans.115 23

[14]Zhang Y,Wang S X,Li S H 1996 Acta Electron.Sin.24 109(in Chinese)[張嚴,王樹勛,李生紅1996電子學報24 109]

[15]Lutes L D,Chen D CK 1991 Int.J.Non-Linear Mech.6 893

[16]Hinich M J1994 Circuits Sys.Signal Proc.13 391

[17]Zhang W J,Huang Y J 2009 Machine Tool and Hydraulics 37 52(in Chinese)[張文俊,黃宜堅2009機床與液壓37 52]

[18]Zheng GB,Jin N D 2009 Acta Phys.Sin.58 4485(in Chinese)[鄭桂波,金寧德2009物理學報58 4485]