田 野，黃仕迪，*

(1. 南方科技大學 力學與航空航天工程系，深圳 518055；2. 南方科技大學 復雜流動及軟物質(zhì)研究中心，深圳 518055)

0 引 言

湍流廣泛存在于各種自然現(xiàn)象和工程應用當中。由于其多尺度耦合、強非線性、強非平衡等復雜特性，關于湍流的研究充滿了挑戰(zhàn)。其中一個焦點是如何通過準確便捷的速度場測量，來同時獲取湍流系統(tǒng)中的大尺度流動結(jié)構和小尺度統(tǒng)計特性，這對于理解湍流的能量級串過程，即能量如何在大、小尺度流動結(jié)構之間傳遞并最終耗散，有著重要的意義[1]。因此，研究人員一直在探索并發(fā)展著各種流體測速技術方法。

目前最為盛行的流體測速技術是粒子成像測速技術（particle image velocimetry，PIV）[2-3]。自20世紀90 年代[4]以來，盡管PIV 在硬件和軟件上都得到了顯著的發(fā)展，但該技術在實際應用中依然存在許多局限，例如無法有效地分辨小尺度流動結(jié)構，也很難準確地計算與強速度梯度相關的物理量等[5-6]。這些局限的根本原因在于，PIV并非直接測量流體速度，而是基于一定大小的問詢域內(nèi)一定數(shù)量的示蹤粒子運動，通過計算相鄰兩張粒子圖像在傅里葉空間的相關性來獲得相應問詢域內(nèi)的流體速度。因此，PIV速度場的準確性和空間分辨率會受限于示蹤粒子濃度尤其是問詢域的大小（通常不小于16 pixel × 16 pixel）。盡管研究人員提出了各種辦法來解決這個困難，包括亞像素插值、多網(wǎng)格迭代、可變形問詢域等，這一由算法原理本身決定的困難依然被認為是PIV的重要挑戰(zhàn)[7-8]。在實際應用中，為了同時獲得大、小尺度流動特性，研究人員一般會采用多臺相機進行全場和局部場的同步拍攝，而這對實驗硬件以及大量數(shù)據(jù)的存儲和后期分析提出了較高的要求。所以，PIV 在同時測量大尺度流動結(jié)構和小尺度統(tǒng)計特性方面并不實用。

另一方面，光流測速法（optical flow velocimetry，OFV）在進行超高空間分辨率的大區(qū)域速度場測量上具有天然的優(yōu)勢[9]。該方法起源于計算機視覺領域，Horn和Schunck[10]基于連續(xù)圖像光強不變的假設，提出了圖像光強亮度限制方程以及相應的變分方程，從而將連續(xù)圖像中剛體運動的速度與光流建立了聯(lián)系。隨后，研究人員采用不同的光流限制方程對OFV進行了完善[11-12]。由于OFV獲得的速度場分辨率與圖像分辨率一致（即單像素分辨率），所以該方法早年也曾被嘗試用于處理PIV的粒子圖像，以期獲得高分辨的速度場[13-15]。然而，不同于剛體運動，流體運動遵循Navier-Stokes方程，而上述的光強限制方程并未考慮流體問題的實際背景，所以光流與流場之間的關聯(lián)缺乏合理的物理依據(jù)。因此，早期的OFV在流體力學領域并未得到關注。

為了解決OFV存在的非物理問題，Liu等[16-18]通過一系列的理論推導和實驗工作，成功地建立了光流與實際流體運動間的定量聯(lián)系，并給出了明確的誤差分析。鑒于OFV在獲取流體速度上有了物理基礎，這個方法逐漸引起了流體力學領域的關注[9]。相比于PIV，OFV可以處理的圖像種類非常豐富，除了粒子圖像，還包括陰影圖和紋影圖[16,19]，甚至是衛(wèi)星云圖[20-21]。需要特別指出的是，陰影圖（紋影圖）等圖像實質(zhì)上是光強變化在投影方向上路徑積分的結(jié)果，所以這類圖像通過OFV得到的結(jié)果一般只能做定性或半定量的分析。雖然粒子圖像可以用OFV進行定量分析，但當粒子濃度不夠高或在相鄰兩幀圖像中的位移過大時，光流方程中的光強時間導數(shù)通常難以準確獲得[16]，從而無法給出準確的流體速度場。

注意到用于PIV的互相關算法對一定范圍內(nèi)的粒子濃度和位移大小并不敏感，Liu和Shen認為可以發(fā)展一種兼容PIV和OFV優(yōu)勢的混合測速法[16]。這個提議近年來被Yang 等[22]和Liu等[23-24]先后獨立實現(xiàn)，在此稱之為基于互相關的光流測速法（crosscorrelation based optical flow velocimetry，cOFV）。簡單地說，cOFV首先利用互相關算法處理粒子圖像，然后將獲得的粗?；俣葓鲎鳛槌跏技s束，最后通過光流法給出高分辨率的速度場（見圖1）?？梢姡琧OFV兼容了PIV的準確性和OFV的高分辨率，具備同時準確測量大、小尺度流動特性的可能性。類似的思路和方法在近年來的其他一些工作亦有體現(xiàn)[25-28]，但大多數(shù)工作都只關注這類方法與傳統(tǒng)方法在獲取流場結(jié)構及速度分布上的比較，很少有研究對其獲取小尺度統(tǒng)計特性的準確性進行仔細檢驗。

圖1 cOFV的分析流程圖Fig. 1 Flow chart of cOFV

本文選取Rayleigh-Bénard（RB）對流作為測試系統(tǒng)，考察cOFV在同時測量大尺度流場結(jié)構和小尺度湍流統(tǒng)計特性上的有效性和準確性。RB對流是通過對系統(tǒng)上、下邊界施加溫差而驅(qū)動流體運動的理想系統(tǒng)。作為一個邊界條件明確的封閉系統(tǒng)，RB對流不僅在實驗和數(shù)值模擬上容易實現(xiàn)和控制，而且兼具壁湍流和各向同性湍流的特性，因此經(jīng)常被用于研究熱對流現(xiàn)象和湍流的基本機理[29-31]。此外，RB對流內(nèi)部的能量耗散率與系統(tǒng)的控制參數(shù)有明確的標度關系[32]，這為進行理論或技術方法的定量檢驗提供了可靠的參考標準。事實上，目前關于RB對流的主流理論[33]的一個重要假設，就是系統(tǒng)內(nèi)部中央?yún)^(qū)域的能量是由大尺度流動結(jié)構注入并最終耗散。為了驗證這個假設，就必須同時測量大、小尺度的流動統(tǒng)計特性。這也是我們選取RB對流作為測試系統(tǒng)考察cOFV的主要動機。

本文往后章節(jié)的結(jié)構如下：我們先介紹本工作中所采用的RB對流實驗裝置、測量技術、以及cOFV的原理和數(shù)據(jù)處理流程；然后展示基于cOFV得到的平均場和流動強度，并與PIV的結(jié)果進行比較；最后展示通過cOFV得到的速度結(jié)構函數(shù)和能量耗散率，并進行討論和總結(jié)。

1 實驗裝置和技術方法

1.1 實驗裝置和測量技術

本研究采用的對流槽與Li等[34]所使用的幾乎相同，是一個長寬高（對應笛卡爾坐標系下的x、y、z三個方向）分別為L= 25.1 cm、D= 7.6 cm、H= 24.8 cm的窄方腔對流槽。由于cOFV對光強變化十分敏感，所以對流槽的導熱銅板表面被均勻地噴涂上了一薄層黑色的鐵氟龍，用以避免表面反光帶來的影響。同時，為了減小周圍環(huán)境的影響，整個對流槽放置在溫度變化小于0.1℃的恒溫箱中。RB 對流的系統(tǒng)控制參數(shù)為衡量浮力驅(qū)動強度的Rayleigh數(shù)（Ra=αgΔTH3/υκ）及反映工作流體物性的Prandtl數(shù)（Pr=υ/κ），其中g是重力加速度，ΔT是系統(tǒng)上、下邊界的溫差，α、υ和κ分別是流體的熱膨脹系數(shù)、運動黏度和熱擴散率。在本實驗中，Ra數(shù)的范圍是9.5×108≤Ra≤4.1×109，而Pr數(shù)保持不變?yōu)?.4。

本實驗獲取粒子圖像的方法與標準的PIV技術一致，在此僅作簡要介紹。我們首先將直徑為50 μm、密度為1.03 g/cm3的示蹤粒子（斯托克斯數(shù)St1×10?4，可以有效地響應流場變化）均勻分布于對流槽中，然后用Nd-YAG激光（波長532 nm）照亮測量區(qū)域（對流槽的整個xz中軸面），再通過Hisense Zyla 相機（選用Nikon 60 mm定焦鏡頭，放大系數(shù)即相機芯片尺寸和圖像尺寸之比為0.049，原始圖像大小為2 000 pixel×2 000 pixel）以4 Hz的雙幀模式采集粒子圖像，其中相鄰兩幀圖像的時間間隔隨著Ra數(shù)的增加由20 ms下降為5 ms。實驗過程中，每個Ra數(shù)下拍攝得到的粒子圖像均大于20 000 張，然后分別采用PIV（Dantec DynamicStudio Ver. 7.0）和cOFV兩種方法進行分析。對于PIV，我們選用的問詢域大小為32pixel×32 pixel，空間重合率為50%，由此獲得的速度矢量間隔為16 pixel×16 pixel，速度矢量個數(shù)為125×125，對應的空間分辨率約為2.1 mm。而對于cOFV，為了避免高Ra數(shù)下粒子位移過大帶來的影響，粒子圖像預先做了降采樣處理至1 000 pixel×1 000 pixel。由于cOFV的空間分辨率為單像素，所以最終獲得的速度矢量個數(shù)為1 000 ×1 000，對應的空間分辨率約為0.27 mm。根據(jù)以往的研究結(jié)果[35-36]和本研究后續(xù)的分析驗證，本實驗的采樣時間小于系統(tǒng)內(nèi)部的Kolmogorov時間（tη1.0 ～ 1.7 s），cOFV的空間分辨率也小于相應的Kolmogorov尺度（η0.91～1.2 mm），滿足進行小尺度統(tǒng)計量分析的要求。

1.2 基于互相關的光流測速法

由于PIV技術已經(jīng)非常成熟，所以在此不作贅述，我們主要介紹cOFV的方法原理。Liu等[17]曾定量給出OFV應用于粒子圖像的誤差來源，主要包括粒子位移、粒子直徑、粒子速度梯度和粒子光照強度。其中，粒子位移與總誤差成正比，且存在一個使總誤差最小的最佳光強。該誤差分析同樣適用于cOFV[24]。據(jù)此，我們在使用cOFV分析粒子圖像之前，首先通過降采樣將初始粒子圖像縮小為原來的一半，從而盡可能降低粒子位移帶來的影響。此外，考慮到相鄰兩幀粒子圖像之間的光強可能存在變化，我們先將二者的平均光強提升（或降低）至相同值，以此控制總體光強不變，再通過高斯濾波使局部光強的變化更為平滑。這些圖像處理可以極大地降低誤差，提高結(jié)果的可信度。

cOFV的分析流程如圖1所示[23]。在完成原始粒子圖像的預處理之后，我們先采用開源程序PIVlab[37]的互相關算法，獲取粗?；俣葓觥；诖舜至；俣葓?，我們通過位移算法（空間插值）獲得粒子圖像1的位移圖像1。鑒于重構后的位移圖像1與粒子圖像2之間的粒子位移已經(jīng)足夠小，我們便可以應用OFV計算二者之間的精細化速度場。為求解光流方程，本文采用Horn-Schunck估計[10]作為初始解，然后借鑒Liu-Shen估計[16,18]計算精細解。最終獲得的準確速度場是空間插值后的粗?；俣葓雠cOFV輸出的精細化速度場的疊加。區(qū)別于Yang等的cOFV[22]，本文的cOFV包含了迭代過程。通常情況下，一次迭代足以獲得準確的速度場[23]。

2 結(jié)果與討論

2.1 平均流場和Reynolds數(shù)

我們首先對比cOFV和PIV獲得的平均速度場（圖2，Ra= 2.0×109）。兩種方法均能清晰地觀察到RB對流的典型流態(tài)[38]：一個呈橢圓形狀的大尺度結(jié)構和兩個反向的角渦。從云圖的顏色分布及色標范圍來看，兩種方法得到的速度分布和大小基本一致，最大速度的差別僅為2%（如圖2各子圖中的色標最大值所示）。鑒于PIV在測量速度場上的應用十分廣泛，而且認可度高，這些半定性半定量的比較在一定程度上說明了cOFV的可靠性。這一結(jié)論和之前的研究一致[22-23,25,27]。

圖2 Ra = 2.0×109下的平均速度場Fig. 2 Time-averaged velocity fields at Ra = 2.0 × 109

為了定量比較cOFV和PIV在獲取平均場以及瞬時場上的區(qū)別，我們考察反映流動強度的Reynolds數(shù)，即Re=VH/υ，其中V為系統(tǒng)內(nèi)部流動的特征速度。對于平均場，選取速度最大值Vm={[(〈U〉t)2+(〈W〉t)2]0.5}max來定義Rem；而對于瞬時場，則選取均方根速度Vrms=[〈(U2+W2)〉V,t]0.5來定義Rerms。其中， 〈···〉t表示時間平均，而 〈 ···〉V,t表示先計算體積平均，再計算時間平均。圖3展示了兩種定義下的Re數(shù)隨著Ra數(shù)的變化情況?？梢钥吹剑琧OFV和PIV獲得的Rem雖然在大小上有微小的系統(tǒng)差別，但它們與Ra數(shù)的標度關系基本一致；而對于Rerms，兩種方法獲得的結(jié)果幾乎完全相同。值得一提的是，本實驗中獲得的Re～Raα標度關系不僅符合關于RB對流的主流理論預測[29,33]，也與以往研究在類似實驗裝置中獲得的結(jié)果基本吻合[39-41]。這些結(jié)果說明了cOFV在定量獲取流場上的準確性，為進行小尺度統(tǒng)計分析提供了基礎。

圖3 Re數(shù)與Ra數(shù)的依賴關系Fig. 3 Log-log plot of Re as a function of Ra

2.2 速度結(jié)構函數(shù)和能量耗散率

為了驗證cOFV獲取小尺度統(tǒng)計量的準確性，我們考察二階速度結(jié)構函數(shù)在湍流耗散區(qū)的行為。前人研究表明[35,42-44]，RB對流中央?yún)^(qū)域的流動統(tǒng)計特性與均勻各向同性湍流相符。對于均勻各向同性湍流，其二階縱向速度結(jié)構函數(shù)DLL和橫向速度結(jié)構函數(shù)DNN在湍流耗散區(qū)滿足以下關系：

其中r為空間兩點的間隔尺度， 〈ε〉為相應區(qū)域內(nèi)的平均能量耗散率[1]。在本研究中，我們選取0.25 ≤x/L≤0.75、0.25 ≤z/H≤ 0.75 的范圍作為中央?yún)^(qū)域，分別計算豎直和水平速度的二階速度結(jié)構函數(shù)，即D[w(r)]=〈[w(x+r)?w(x)]2〉和D[u(r)]=〈[u(x+r)?u(x)]2〉，其 中w=W?〈W〉t和u=U?〈U〉t分別為豎直和水平方向上的脈動速度。

我們首先比較基于cOFV和PIV獲得的二階速度結(jié)構函數(shù)（Ra= 2.0×109）。圖4展示了根據(jù)公式（1）進行數(shù)據(jù)補償后的結(jié)果，其中橫坐標采用Kolmogorov尺度ηc=υ0.75/( 〈ε〉c)0.25進行了無量綱化，所用的 〈ε〉c詳見后文平均能量耗散率的比較分析?？梢钥吹?，當r/ηc＞ 10時， cOFV和PIV獲得的結(jié)果幾乎一致；當2 ＜r/ηc＜ 10時，PIV給出的結(jié)果明顯偏大，這和PIV實驗測量在小尺度上的噪聲有關（見參考文獻[30]中的圖7和相關討論）；而當r/ηc＜ 2時，PIV已經(jīng)無法提供數(shù)據(jù)，cOFV的結(jié)果則按理論預測的趨勢逐漸形成一個“平臺”?；谶@些比較可知：一方面，PIV的空間分辨率不足以解析湍流耗散尺度；另一方面，cOFV不僅能解析耗散區(qū)的統(tǒng)計特性，而且降采樣、高斯濾波等圖像預處理技術有效地降低了實驗測量的噪聲，更適合進行小尺度統(tǒng)計量的分析。因此下文我們僅討論cOFV的結(jié)果。

圖4 Ra = 2.0×109下二階速度結(jié)構函數(shù)Fig. 4 Second-order velocity structure functions at Ra = 2.0×109

圖5和圖6分別展示了cOFV獲得的不同Ra數(shù)下豎直和水平速度的二階速度結(jié)構函數(shù)。這些結(jié)果也按圖4的方式進行了補償和無量綱化處理。從結(jié)構函數(shù)的形態(tài)上看，D[w(r)] 和D[u(r)]在Kolmogorov尺度以下都存在一定的“平臺”，但不同于D[w(r)]的縱向和橫向分量在湍流耗散區(qū)幾乎重合，D[u(r)]的兩個分量有一定程度的錯開。這說明在本實驗裝置中，中央?yún)^(qū)域的豎直脈動速度更接近均勻各向同性湍流的特性，而水平脈動速度則有可能受到了窄方腔對流槽的空間約束的影響[45]。盡管如此，我們依然可以結(jié)合公式（1），通過耗散區(qū)的“平臺”（如圖5和圖6中的水平實線所示，其中圖6的水平實線是相同Ra數(shù)下縱向和橫向結(jié)構函數(shù)對應“平臺”的平均）獲得相應的平均能量耗散率 〈ε〉SFw,c和 〈 ε〉SFu,c。

圖5 不同Ra數(shù)下豎直速度二階結(jié)構函數(shù)與r的依賴關系Fig. 5 Log-log plot of second-order structure functions of vertical velocity versus r for different Ra

圖6 不同Ra數(shù)下水平速度二階結(jié)構函數(shù)與r的依賴關系Fig. 6 Log-log plot of second-order structure functions of horizontal velocity versus r for different Ra

另一方面，我們也可以根據(jù)定義直接計算能量耗散率，并和基于速度結(jié)構函數(shù)獲得的能量耗散率進行對比驗證。考慮到本實驗中獲得的速度場是一個二維場，所以我們需要通過現(xiàn)有的二維數(shù)據(jù)重構出三維的耗散場。為此，我們借鑒Doron等[46]和Luznik等[47]的做法，假設三維耗散場中缺失的項與二維耗散場中對應的項相近，由此得到的能量耗散率重構公式如下：

圖7展示了通過公式（2）重構的能量耗散率場（Ra= 2.0×109）。根據(jù)疊加的平均速度場，我們可以看到中央?yún)^(qū)域的能量耗散率明顯小于靠近邊界區(qū)域的強度，且耗散率較大的地方主要集中在大尺度流動結(jié)構即剪切較強的區(qū)域。需要指出的是，由于RB系統(tǒng)在邊界附近存在強烈的溫度漲落，其引起的局部光強不均勻使得相應區(qū)域的光折射率變化極大，進而導致局部光強非常不均勻。而cOFV對光強的變化十分敏感，在計算邊界附近的速度場尤其是速度梯度時會存在較大誤差。為了避免邊界效應帶來的影響，下文只對RB系統(tǒng)中央?yún)^(qū)域（即遠離邊界處）的能量耗散率進行定量分析。

圖7 Ra = 2.0×109下的能量耗散率場Fig. 7 The reconstructed energy dissipation field at Ra = 2.0×109

根據(jù)重構后的能量耗散場，我們便可以得到中央?yún)^(qū)域的平均能量耗散率 〈ε〉D,c，并與基于速度結(jié)構函數(shù)間接獲得的平均能量耗散率 〈ε〉SFw,c和 〈ε〉SFu,c進行比較。從圖8可以看到，三種平均能量耗散率在誤差范圍內(nèi)幾乎完全相同，說明它們在有效性上可以互為支撐驗證。

為了檢驗這些實驗結(jié)果的正確性，我們選取對流槽幾何參數(shù)與本研究接近的直接數(shù)值模擬的結(jié)果[44]一起畫在圖8中作為參考標準。需要指出的是，對于RB對流這類封閉系統(tǒng)，系統(tǒng)的能量耗散率和傳熱效率有明確的對應關系，而能量耗散率的空間分布則和大尺度流動結(jié)構有關。以往的研究[44-45,48-49]表明，RB系統(tǒng)在窄方腔和立方體對流槽中的傳熱效率基本相同，大尺度流動結(jié)構也非常相似，并沒有表現(xiàn)出明顯的空間約束效應。因此，采用立方體RB系統(tǒng)的數(shù)值模擬結(jié)果進行對比，具有合理性?？梢钥吹?，本研究的實驗結(jié)果和數(shù)值結(jié)果符合得很好，除了Ra= 9.5×108那組數(shù)據(jù)有所偏離?？赡茉蚴菍τ谶@個Ra數(shù)，統(tǒng)一選取的中央?yún)^(qū)域（0.25 ≤x/L≤0.75，0.25 ≤z/H≤ 0.75）偏大，導致部分大尺度流動結(jié)構亦有貢獻。如果對圖8中的結(jié)果進行擬合，我們得到的標 度 關 系 為 〈ε〉c～Ra1.31±0.04?；?顧 引 言 所 述，關于RB對流的主流理論[33]假設，系統(tǒng)中央?yún)^(qū)域的能量耗散率與大尺度流動結(jié)構的能量注入相平衡，即〈ε〉GL,c～V3/H=υ3Re3/H4。鑒于此處考慮的是空間平均量，我們代入前文中獲得的標度關系Rerms～Ra0.46，由此得到的理論預測為 〈ε〉GL,c～Ra1.38。這與圖8的擬合結(jié)果非常接近，一方面支持了理論假設，另一方面也表明通過cOFV同時獲得大、小尺度流動特性具有可靠性。

圖8 中央?yún)^(qū)域的平均能量耗散率與Ra數(shù)的依賴關系Fig. 8 The energy dissipation rate in the central region as a function of Ra

3 結(jié) 論

本文采用基于互相關的光流測速法（cOFV），對窄方腔RB對流系統(tǒng)中的大尺度流動結(jié)構和小尺度湍流統(tǒng)計特性進行了同時測量和研究。實驗中的Pr數(shù)保持不變?yōu)?.4，Ra數(shù)的變化范圍是9.5×108≤Ra≤ 4.1×109，主要結(jié)果如下：

1）在大尺度流動結(jié)構方面，cOFV與PIV給出的流場形態(tài)相同。定量地，兩種測速方法給出的流場強度互相吻合，其中基于均方根速度定義的Rerms數(shù)無論在數(shù)值大小還是與Ra數(shù)的標度關系上都完全一致，相應的標度指數(shù)也與前人的實驗和理論工作相符合。這些結(jié)果表明，cOFV在測量全局場上具備與PIV技術同等的能力。

2）為了探究cOFV在測量小尺度湍流統(tǒng)計量上的準確性，我們考察了RB系統(tǒng)中央?yún)^(qū)域的速度結(jié)構函數(shù)和能量耗散率。我們根據(jù)均勻各向同性湍流在耗散區(qū)的特性，通過二階速度結(jié)構函數(shù)間接得到了能量耗散率，并與根據(jù)定義直接計算獲得的能量耗散率進行對比，發(fā)現(xiàn)二者和已有的直接數(shù)值模擬結(jié)果相吻合，且它們與Ra數(shù)的標度關系也和主流理論預測一致。

這些結(jié)果表明基于互相關的光流測速法在同時測量大、小尺度流動特性上有很好的應用潛力。我們今后將會把此方法應用到其它湍流系統(tǒng)中，通過把不同尺度的流動結(jié)構和能量耗散率的空間分布建立關聯(lián)，來檢驗湍流理論的一些重要假設[1]。

致謝：感謝劉天舒教授把基于互相關的光流測速法介紹給我們，并提供了開源程序（https://github.com/Tianshu-Liu）；同時也感謝高振源對實驗測量的幫助，以及黃永祥教授的相關討論。