兩相流有限穿越可視圖演化動力學研究*

高忠科 胡瀝丹 周婷婷 金寧德

(天津大學電氣與自動化工程學院，天津 300072)

(2012年12月27日收到;2013年1月27日收到修改稿)

1 引言

在石油、化工及核反應堆等工業領域存在著大量的兩相流動現象，例如動力工程中鍋爐蒸發管中的蒸汽-水兩相流動，化工工業中物料輸送管道和反應釜攪拌器中的氣液兩相流以及石油工程中油田開采油井中的油水兩相流和氣液兩相流等.在兩相流動過程中，由于相界面的形狀和分布隨著時間和空間的改變而變化，且相間存在不可忽略的相對速度，導致相界面呈現出不同的幾何形狀或流動結構，稱之為兩相流流型.流型不僅影響混合流體的流動特性、壓降特性和傳熱傳質性能，而且還是兩相流測量、控制與預測系統的重要參數.因此，揭示兩相流不同流型的內在流動特性對兩相流工業系統優化設計及工況動態監測具有重要實際意義.兩相流系統是具有混沌、耗散、有序與無序等復雜特征的非線性動力學復雜系統[1-3]，由于兩相流受到諸如流體湍動、相間界面相互作用及相間局部相對運動等許多復雜因素的影響，其流動行為呈現出高度無規則性、多態性和結構不穩定性，至今其流型轉變動力學機理尚未十分清楚.兩相流流型的實驗測量、動態描述、結構識別及演化動力學機理揭示等相關研究已成為多相流研究領域的一個重要方向.

本文針對小管徑垂直上升管氣液兩相流，全新優化設計了弧形對壁式電導傳感器，通過動態實驗在獲取傳感器測量信號的基礎上，采用我們最近提出的有限穿越可視圖理論分析測量信號揭示不同氣液兩相流流型的內在演化動力學機理.可視圖理論屬于復雜網絡時間序列分析理論的一種，復雜網絡理論是對復雜系統的一種抽象和描述方式，所謂網絡就是由節點和連邊組成，其中節點表示系統的元素，兩節點的連邊表示元素之間的相互作用.在小世界網絡[4]和無標度網絡[5]的開創性工作之后，復雜網絡理論已成功應用于物理學、生物學、經濟學和工程學等多個學科領域的大量復雜系統研究之中，并在相應領域取得了很大進展.近年來，基于時間序列測量信號的復雜網絡動力學研究受到了不同領域研究人員的廣泛關注[6-19]，不同的時間序列復雜網絡構建與分析算法相繼提出[20-33]并已成功應用于不同的復雜動力學系統研究之中[34-48].值得指出的是，Lacasa和Luque[24-26]等提出的可視圖和水平可視圖建網方法具有可處理數據量大且計算速度快的優點，從而具備潛在的實時測量數據在線分析能力.我們最近研究發現這兩種可視圖建網方法雖然可以有效區分和刻畫不同類型的動力學信號，但其在實測信號分析中的噪聲適應性上存在問題.為此，在可視圖建網思想啟發下，我們提出了一種有限穿越可視圖建網方法[32]，通過對三種可視圖建網方法(可視圖、水平可視圖、有限穿越可視圖)的比較評價，發現有限穿越可視圖不僅可以有效區分具有不同動力學特性的時間序列，同時在兩相流實測信號分析上表現出較好的抗噪能力.因此，本文采用有限穿越可視圖理論研究小管徑氣液兩相流流型演化動力學機理.此外，需要指出的是我們前期針對內徑為125mm的大管徑兩相流流動特性，提出了不同類型的兩相流復雜網絡，并將其應用于兩相流流型非線性動力學特性研究之中[2,41,42].與前期研究不同，與大管徑(125 mm直徑)兩相流相比小管徑(20 mm直徑)氣液兩相流流型動力學演化機理更為復雜，本文針對小管徑垂直上升氣液兩相流全新設計弧形對壁式電導傳感器并通過動態實驗獲取傳感器測量信號，在此基礎上采用具有可處理數據量大、計算速度快、抗噪能力強等優點的有限穿越可視圖理論分析測量信號，研究不同氣液兩相流流型的非線性動力學特性，尤其是在不同流型演化過程中的細節演化動力學機理.用傳統的科學觀念和方法研究復雜的兩相流動問題仍存在較大局限性，本文從全新理論產生及信息處理技術不斷完善的角度采用有限穿越可視圖理論認識兩相流動現象，從而為研究復雜混相流體流動結構及流動特性提供了一種新的理論工具.

2 弧形對壁式電導傳感器設計與兩相流動態實驗數據獲取

為了有效獲取內徑為20 mm的垂直細管氣液兩相流不同流型的流動信息，本文采用有限元分析軟件ANSYS優化設計了弧形對壁式電導傳感器，傳感器實物結構圖如圖1所示，其材質為不銹鋼，其中圖示半環形電極成對內嵌在測量管道內，構成一組檢測電極和兩組保護電極，每組電極處于管道徑向的同一個截面上且嚴格對稱，管道外側添加保護層增加傳感器機械強度.我們在天津大學油氣水多相流流動實驗裝置上進行了垂直氣液兩相流的動態實驗.油氣水多相流流動實驗裝置可以實現在20 mm管徑的樹脂管道中進行氣相、水相、油相及聚合物的混合流動實驗，控制不同相的流量，并對其進行測量.動態實驗測量系統主要包括信號調理模塊、數據采集模塊和數據存儲模塊.信號調理模塊主要由C/V轉換電路，整流電路，濾波電路和放大電路四部分構成.數據采集模塊選用美國國家儀器(NI)公司基于PXI總線技術的PXI 4472和PXI 6221數據采集卡進行數據采集，采樣頻率為4000 Hz，采樣時間為30 s.數據存儲模塊采用PC機，通過數據采集卡匹配的圖形化編程軟件LabVIEW實現，可以對數據進行實時的波形顯示和實時存儲等功能.

圖1 弧形對壁式電導傳感器實物圖

實驗方案為首先固定一個水相流量，然后逐漸增加氣量，當氣液兩相流量均穩定后，觀察并記錄該工況的流型并利用弧形對壁式電導傳感器測量得到每一個工況下的電壓波動信號.改變水相流量，重復上述過程，得到一系列對應于不同流型的時間序列測量信號.實驗過程中，隨著工況的改變，兩相流流型也會隨之變化.本次實驗一共采集到對應于泡狀流、段塞流和混狀流三種典型流型的多組實驗數據，其中水相流量變化范圍為0.34—1.29 m/s，氣相閥門開度的變化范圍為5—25 psi(1 psi=6.895 kPa)，水相密度為1.0 g/cm-3.在正弦信號的激勵下，氣相的變化引起檢測電極上的電壓信號波動.對應于不同氣液兩相流流型的弧形對壁式電導傳感器測量信號如圖2所示.圖2(a)為細小泡狀流，泡徑約1 mm，此時測量信號的細小波動反映了該流型下的小氣泡隨機運動特性;圖2(b)為泡狀流，泡徑約2—3 mm，此時隨著小氣泡泡徑的增大，測量信號波動變強，表明該傳感器能夠反映出因泡徑變化所引起的流型流動特性的變化;圖2(c)為段塞流，液塞與氣塞交替出現，高電平信號表示液塞通過測量截面，高電平信號的微小波動表明水相中夾帶小氣泡，負脈沖信號表明有氣塞通過測量截面;增加氣相流量，出現圖2(d)所示的段塞流，其負脈沖變長，表明隨著氣相流量增大，氣塞變長;圖2(e)為混狀流，此時測量信號上下振動劇烈，表明該流型中氣相與液相上下震蕩的不穩定運動特性;增加氣相流量，出現圖2(f)所示的混狀流，此時測量信號以低電平為主，表明在此流動工況下氣相運動占優勢.因此，弧形對壁式電導傳感器對氣液兩相流型有較高的空間分辨率，能夠刻畫不同流型下的流體流動細節信息，為深入研究氣液兩相流流型演化動力學機理提供了可靠有效的信息源.

圖2 弧形對壁式電導傳感器測量信號 (a)泡狀流;(b)泡狀流;(c)段塞流;(d)段塞流;(e)混狀流;(f)混狀流

3 有限穿越可視圖建網分析算法

受可視圖建網方法思想的啟發，我們最近提出一種有限穿越可視圖建網分析算法，通過對可視圖、水平可視圖和有限穿越可視圖的比較評價，發現有限穿越可視圖不僅可以有效區分具有不同動力學特性的時間序列，同時在兩相流實測信號分析上表現出較好的抗噪能力.我們如下將簡要介紹可視圖和有限穿越可視圖建網分析算法.

3.1 可視圖建網方法

Lacasa等[24]提出了可視圖建網方法(Visibility Graph).對于離散時間序列，將數據點定義為網絡節點，數據點之間滿足可視性準則的連線定義為網絡連邊.如圖3所示，用直方條表示一個周期時間序列中20個數據點(圖3(a))，數據真實值與直方條高度相對應.若兩個直方條頂端相互可視，則認為兩點在網絡中相連.圖3(b)是由該方法生成的網絡，依次排開的實點與離散時間序列數據點一一對應，實點之間的連線與直方條之間的可視線一一對應.值得指出的是，可視線不能重復或者遺漏，節點與自身不能相連，兩個直方條之間的可視線不能穿越其他直方條.

可視性準則：如果離散時間序列中點(ta,ya)和點(tb,yb)相連接，那么對于任意的點(tc,yc)，其中ta＜tc＜tb，滿足

可視圖網絡具有如下性質：每一個點至少和它的左鄰點和右鄰點相連;網絡為一個無向無權網絡;橫軸和縱軸坐標尺度變化或者經過仿射變換后，可視性保持不變.可視圖方法能夠繼承原始時間序列的一些特征，即周期時間序列轉化成規則網絡，隨機時間序列轉化成隨機網絡，分形時間序列轉化成無標度網絡.

圖3 時間序列可視圖建網示意圖

3.2 有限穿越可視圖建網方法

在可視圖建網方法思想的啟發下，我們最近提出了一種有限穿越可視圖建網方法[32](Limited Penetrable Visibility Graph).其基本思想是：定義有限穿越視距N，基于可視圖方法可視性準則，規定若網絡中的兩節點相連，當且僅當與節點相對應的兩個直方條之間的連線被位于其間的其他直方條截斷次數n≤N.圖4為一個有限穿越可視圖建網示例，有限穿越視距N=1，圖4(a)中，實連接線表示基于可視圖方法的連接線，虛連接線表示新增連接線，圖4(b)中，實點上方為基于可視圖方法的連接線，下方為新增連接線.

有限穿越可視圖方法具有連接無向性與仿射變換或者坐標軸尺度變換后可視性不變的特點，每個節點至少和2(N+1)個節點相連.研究發現，在噪聲的干擾下，可視圖方法中本應相連的節點易受噪聲點隔斷，對建網效果產生影響，動力學特性被噪聲掩蓋.由圖3和圖4可知，有限穿越可視圖網絡具有更強的連接性，長程連接能夠反映序列波動趨勢，在有噪聲的情況下仍能較好的反映信號的動力學特性.

4 氣液兩相流有限穿越可視圖流型動力學分析

實驗中采用弧形對壁式電導傳感器測取了對應于不同流動工況的23組測量信號，其中10組測量信號對應于泡狀流，4組測量信號對應于段塞流，9組測量信號對應于混狀流.選取有限穿越視距N=5，對實驗測取的23組電壓波動信號進行有限穿越可視圖建網分析，獲得網絡度分布曲線.定義網絡節點平均度值為avK，曲線峰值記為Nmax(即具有相同度值節點個數的最大值)，以avK為橫坐標，Nmax為縱坐標，得到如圖5所示的氣液兩相流有限穿越可視圖網絡特征參數分布圖.由圖可見，泡狀流平均度值較小，且集中分布在80—130之間，Nmax分布較廣泛，在40—180之間;混狀流平均度值較高，且分布較廣泛，在200—450之間，Nmax較小，在10—60之間;段塞流平均度值和Nmax主要分布在泡狀流和混狀流的中間區域.我們選取三種流型中有代表性的點進一步揭示不同流型的網絡度分布狀態.泡狀流中選取Nmax較高的工況點(a)和Nmax較低的工況點(b)，段塞流選取工況點(c)，混狀流選取平均度值較低的工況點(d)和平均度值較高的工況點(e)，如圖5所示.圖6為上述所選的對應于不同流型5個工況點的電壓波動信號和相應的有限穿越可視圖度分布曲線.

圖4 有限穿越可視圖建網示意圖

圖5 氣液兩相流流型點在avK-Nmax特征參數平面分布

圖6(a)流型為泡狀流，其中氣泡泡徑約為2—3 mm，測量信號以高電平為主且具有明顯的波動.圖6(b)流型為細小泡狀流，氣泡泡徑約為1 mm，測量信號有細小的波動.通過比較分析泡狀流流型特征和Nmax值發現，隨著Nmax值的減小，即由流動工況(a)到流動工況(b)方向，泡徑呈現減小趨勢，表明有限穿越可視圖網絡指標可有效刻畫兩相流流動物理結構細節特征.圖6(c)流型為段塞流，其電壓波動信號具有擬周期特性，相應度分布圖由多個尖峰構成.雖然(b)，(c)兩點在avK-Nmax特征參數平面位置接近(如圖5所示)，但是其度分布圖截然不同(如圖6(b)，(c)所示)，表明不同流型具有明顯不同的度分布特征.圖6(d)流型為混狀流，其電壓波動信號具有明顯的上下跳躍特征，因此節點之間相互“遮擋”，節點可視性相對較低.圖6(e)流型為減小水相流速時的混狀流，此時電壓波動信號震蕩頻率變緩，節點之間“遮擋”減弱，節點具有相對較高的可視性，相應其avK值變大.通過比較分析混狀流流型特征和avK發現，隨著avK值的增大，即由流動工況(d)到流動工況(e)方向(如圖5所示)，水相流速呈增大趨勢，表明有限穿越可視圖網絡指標可有效表征兩相流流動物理結構的宏觀特征.

為了進一步揭示兩相流流型演化動力學特性，我們從有限穿越可視圖中計算對應于不同流型的異速生長指數，以期通過異速生長指數刻畫流型演化動力學復雜性.異速生長尺度的概念最早是在生物學研究領域由West等人于1997年和1999年提出的[49,50]，異速生長簡單的說就是指不同的身體部分存在著不同的生長速率，尺度關系則是生態學中的核心概念，它以簡單的形式量化復雜的系統的一般性特征.在生物學研究中，所有的生物尺度可以用如下異速生長方程來描述：

式中Y為表征有機體生理特征的一個變量，例如新陳代謝速率等，M則是指有機體的大小或者體重，b為尺度冪指數.在許多實際復雜網絡研究中都能夠發現異速生長尺度率，如食物網絡[51]，器官代謝網絡[52]，世界貿易投資網絡[53,54]等.復雜網絡中的異速生長尺度率主要依賴于網絡的二叉樹結構.一棵樹擁有著一個根節點和許多的分支以及葉節點，所以我們可以把樹看作是從根節點到葉節點的所用分支和節點的集合.對于樹中的每一個節點i，我們將它的Ai值定義為i節點下屬的所有葉節點A值的和再加1，顯然最底層葉節點沒有子節點，所以它們的A值統一定義為1，計算各節點的A值時是從樹的最底層葉節點開始依次遞歸計算各節點A值，直到根節點為止，如下式所示：

其中 j是i的子節點.樹中的每一個節點的C值定義與A值定義相似，我們定義Ci為i節點下屬的所有葉節點C值的和然后再加上i節點的A值，其中最底層的節點C值也統一定義為1，計算各節點的C值的方法與A值相似，也是從下往上逐層遞歸計算，迭代計算過程如下式所示：

其中 j是i的子節點.當一股流從樹的根節點到葉節點，那么A值代表著節點的尺寸規模，C值代表著流向節點的流量.根據異速生長方程可知，如果該樹滿足異速生長尺度關系，那么流量C和尺寸A應該滿足冪律關系，即有

其中冪率指數η是網絡樹的一個很重要的結構參數.鏈狀樹和星型樹是兩種極端的樹，通常冪率指數η會介于二者之間，對于鏈狀樹其指數η=2，對于星型樹其指數η=1，所以異速生長冪率指數范圍為：1＜η＜2.研究發現[55]，不是所有的樹都滿足異速生長尺度關系，對于同一個網絡，構造的分岔樹不同那么對應的異速生長尺度關系也會不同.一個復雜網絡可以構造出不同的分岔樹結構，本文從復雜網絡中構造的是最小分岔樹結構，具體構造過程為：首先通過有限穿越可視圖算法得到與時間序列對應的復雜網絡，然后找出該網絡中度值最大的點作為最小分岔樹的根節點，接著找出所有與根節點相連的點作為它的子節點，依次類推，直到整個復雜網絡完全轉化為最小分岔樹為止.通過計算分析，我們發現對應于不同氣液兩相流流型的有限穿越可視圖復雜網絡最小分岔樹結構都滿足異速生長尺度關系，即AC曲線在雙對數坐標下滿足冪律分布的形式，而可視圖復雜網絡對應的最小分岔樹結構卻不具有很好的異速生長尺度關系，即AC曲線在雙對數坐標下不具有單一的線性關系，如圖7所示，表明有限穿越可視圖在氣液兩相流異速生長指數分析上的優越性.

圖6 典型工況點的電壓波動信號和有限穿越可視圖度分布曲線

圖7 不同流型的有限穿越可視圖及可視圖網絡AC數值分布圖 (a)，(c)，(e)為有限穿越可視圖對應于三種典型流型的AC數值分布圖;(b)，(d)，(f)為可視圖對應于三種典型流型的AC數值分布圖

圖8 不同流型的有限穿越可視圖網絡異速生長指數和網絡平均度值的聯合分布圖

對測量的23組電壓波動信號進行有限穿越可視圖建網并提取相應異速生長指數，異速生長指數η和網絡平均度值avK的聯合分布如圖8所示.由圖可見，泡狀流、段塞流和混狀流流型點分布在avK-η指標參數平面的不同區域，表明有限穿越可視圖網絡異速生長指數與度分布參數組合對三種流型具有較好識別能力.我們進一步采用自適應最優核(AOK)分布算法[56,57]考察測量信號的時頻分布特征，對應于不同流型的AOK時頻分布圖如圖9所示.泡狀流中存在著大量的小氣泡，泡群運動隨機可變，其動力學行為十分復雜，使得泡狀流運動頻率分布較寬(主要分布在0—40 Hz)，但運動能量較低(集中在0—12范圍內)，如圖9(a)所示;在流型從泡狀流逐漸向段塞流轉變過程中，小氣泡的密度不斷增大，泡群開始發生聚并形成氣塞，氣塞與液塞的擬周期交替運動使段塞流動力學行為變得相對簡單，如圖9(b)所示，與泡狀流相比，段塞流頻率分布較窄(主要集中于0—10 Hz)，能量較高(分布在0—100范圍內).相應由圖8可見，流型在由泡狀流向段塞流演化過程中，網絡異速生長指數逐漸增大，在段塞流時達到最大值，很好地表征了在泡狀流到段塞流演化過程中的動力學特性，網絡異速生長指數越小反映了流體運動頻率范圍越寬、運動能量越低、流型的動力學行為越復雜.隨著氣相流量的進一步增加，流型從段塞流轉變為混狀流，由于流體湍動能量的增強，液體被氣體沖碎形成液塊，受重力作用在整體上升的同時上下翻滾，使得混狀流動力學復雜性高于段塞流，如圖9(c)所示，混塞流運動頻率高于段塞流，其主要分布在0—20 Hz，運動能量低于段塞流，主要集中于0—70范圍內.相應流型在由段塞流向混狀流演化過程中，隨著網絡異速生長指數逐漸減小，表征了氣液兩相流在由段塞流向混狀流演化過程中流體動力學行為復雜性逐漸增大.

圖9 典型流型的電導波動信號自適應最優核時頻分布圖 (a)泡狀流;(b)段塞流;(c)混狀流

因此，通過上述分析可知，有限穿越可視圖的網絡異速生長指數和網絡平均度值的聯合分布可實現對小管徑氣液兩相流的流型辨識;有限穿越可視圖的度分布曲線峰值可有效刻畫兩相流流動物理結構的細節特征;網絡平均度值可表征流動結構的宏觀特性;網絡異速生長指數可揭示兩相流流型演化動力學復雜性，為深入地理解氣液兩相流相間相互作用機理及流動特性提供了參考.

5 結論

針對小管徑垂直上升氣液兩相流全新設計弧形對壁式電導傳感器，通過動態實驗獲取對應于不同流型的傳感器測量信號.采用我們最近提出的有限穿越可視圖理論從測量信號中構建復雜網絡，通過分析構建的對應于不同流動工況的有限穿越可視圖網絡，發現網絡度分布曲線峰值Nmax對兩相流泡徑大小分布十分敏感，可有效刻畫兩相流流動物理結構細節特征;網絡平均度值avK隨著水相流速的增大而逐漸增大，可有效表征兩相流流動物理結構的宏觀特性;網絡異速生長指數η對流體動力學復雜性十分敏感，可揭示不同流型演化過程中的細節演化動力學特性;網絡異速生長指數和網絡平均度值的聯合分布可實現對小管徑氣液兩相流的流型辨識.兩相流測量信號的有限穿越可視圖分析可為兩相流相間相互作用機理揭示及流動特性分析提供新視角、建立新途徑.

[1]Daw C S，Finney C E A，Vasudevan M，van Goor N A，Nguyen K，Bruns D D，Kostelich E J，Grebogi C，Ott E，Yorke J A 1995 Phys.Rev.Lett.75 2308

[2]Gao Z K，Jin N D 2009 Phy.Rev.E 79 066303

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

[4]Watts D J，Strogatz S H 1998 Nature 393 440

[5]Strogatz S H 2001 Nature 410 268

[6]Albert R，Barab′asi A L 2002 Rev.Mod.Phys.74 47

[7]Wang W X，Wang B H，Hu B，Yan G，Ou Q 2005 Phys.Rev.Lett.94 188702

[8]Zheng D F，Hui P M，Trimper S，Zheng B 2005 Physica A 352 659

[9]Pan Z F，Li X，Wang X F 2006 Phys.Rev.E 73 056109

[10]Huang L，Park K，Lai Y C，Yang L，Yang K Q 2006 Phys.Rev.Lett.97 164101

[11]Li X G，Gao Z Y，Li K P，Zhao X M 2007 Phys.Rev.E 76 016110

[12]Pei W D，Chen Z Q，Yuan Z Z 2008 Chin.Phys.B 17 0373

[13]Zhou J，Lu J A，Lu J H 2008 Automatica 44 996

[14]Lu X B，Wang X F，Fang J Q 2010 Physica D 239 341

[15]Zhou T，Kuscsik Z，Liu J G，Medo M，Wakeling J R，Zhang Y C 2010 Proc.Natl.Acad.Sci.USA 107 4511

[16]Guo J L 2010 Chin.Phys.B 19 120503

[17]Duan Z S，Chen G R 2012 Chin.Phys.B 21 080506

[18]Zhang C，Sheng H Z，Li F，Yang H Q，Feng G L 2012 Acta Phys.Sin.61 148902(in Chinese)[張聰，沈惠璋，李峰，楊何群2012物理學報61 148902]

[19]Wang W X，Yang R，Lai Y C，Kovanis V，Grebogi C 2011 Phys.Rev.Lett.106 154101

[20]Zhang J，Small M 2006 Phys.Rev.Lett.96 238701

[21]Zhou C S，Zemanov?L，Zamora G，Hilgetag C C，Kurths J 2006 Phys.Rev.Lett.97 238103

[22]Zhang J，Sun J F，Luo S D，Zhang K，Nakamura T，Small M 2008 Physica D 237 2856

[23]Yang Y，Yang H J 2008 Physica A 387 1381

[24]Lacasa L，Luque B，Ballesteros F，Luque J，Nuno J C 2008 Proc.Natl.Acad.Sci.USA 105 4972

[25]Luque B，Lacasa L，Ballesteros F，Luque J 2009 Phys.Rev.E 80 046103

[26]Lacasa L，Toral R 2010 Phys.Rev.E 82 036120

[27]Xu X，Zhang J，Small M 2008 Proc.Natl.Acad.Sci.USA 105 19601

[28]DongesJF，HeitzigJ，DonnerRV，KurthsJ2012Phys.Rev.Etextbf85 046105

[29]Zou Y，Donner R V，Wickramasinghe M，Kiss I Z，Small M，Kurths J 2012 Chaos 22 033130

[30]Gao Z K，Jin N D 2009 Chaos 19 033137

[31]Gao Z K，Jin N D 2012 Nonlinear Analysis：Real World Applications 13 947

[32]Zhou T T，Jin N D，Gao Z K，Luo Y B 2012 Acta Phys.Sin.61 030506(in Chinese)[周婷婷，金寧德，高忠科，羅躍斌2012物理學報61 030506]

[33]Shimada Y，Ikeguchi T，Shigehara T 2012 Phys.Rev.Lett.109 158701

[34]Fang X L，Jiang Z L 2007 Acta Phys.Sin.56 7330(in Chinese)[方小玲，姜宗來2007物理學報56 7330]

[35]Dong Z，Li X 2009 Acta Phys.Sin.59 1600(in Chinese)[董昭，李翔2009物理學報59 1600]

[36]Li X，Dong Z 2011 Phys.Rev.E 84 062901

[37]Zhou L，Gong Z Q，Zhi R，Feng G L 2008 Acta Phys.Sin.57 7380(in Chinese)[周磊，龔志強，支蓉，封國林2008物理學報57 7380]

[38]Gong Z Q，Wang X J，Zhi R，Feng A X 2011 Chin.Phys.B 20 079201

[39]Gong Z Q，Zhi R，Hou W，Wang X J，Feng G L 2012 Acta Phys.Sin.61 029202(in Chinese)[龔志強，支蓉，侯威，王曉娟，封國林2012物理學報61 029202]

[40]Gao X Y，An H Z，Fang W 2012 Acta Phys.Sin.61 098902(in Chinese)[高湘昀，安海忠，方偉2012物理學報61 098902]

[41]Gao Z K，Jin N D 2009 Chin.Phys.B 18 5249

[42]Gao Z K，Jin N D，Wang W X，Lai Y C 2010 Phy.Rev.E 82 016210

[43]Zhang J，Zhang K，Feng J F，Small M 2010 PLoS Comput.Biol.6 e1001033

[44]Liu C，Zhou W X，Yuan W K 2010 Physica A 389 2675

[45]Shao Z G 2010 Appl.Phys.Lett.96 073703

[46]Gao Z K，Jin N D 2011 Chem.Eng.Sci.66 2660

[47]Walker D M，Tordesillas A 2012 Phys.Rev.E 85 011304

[48]Ghaffari H O Young R P 2012 Europhys.Lett.98 48003

[49]West G B，Brown J H，Enquist B J 1997 Science 276 122

[50]West G B，Brown J H，Enquist B J 1999 Science 284 1677

[51]Garlaschelli D，Caldarelli G，Pietronero L 2003 Nature 423 165

[52]Santillan M 2003 Eur.Phys.J.B 223 249

[53]Duan W Q 2007 Eur.Phys.J.B 59 271

[54]Song D M，Jiang Z Q，Zhou W X 2009 Physica A 388 2450

[55]Ni X H，Jiang Z Q，Zhou W X 2009 Phys.Lett.A 373 3822

[56]Jones D L，Baraniuk R G 1995 IEEE Transactions on Signal Processing 43 2361

[57]Du M，Jin N D，Gao Z K，Wang Z Y，Zhai L S 2012 International Journal of Multiphase Flow 41 91