基于實驗數據的湍流擴散POD模態分析

王芳，賈勝坤，張會書，袁希鋼，余國琮

（化學工程聯合國家重點實驗室，天津大學化工學院，天津 300072）

引言

物質在湍流條件下的擴散是廣泛存在于化工過程中的基本現象[1-2]，了解湍流中的傳質擴散規律，對于化工基礎和應用研究具有重要意義。然而由于其復雜性，湍流條件下物質擴散或濃度分布的嚴格模擬十分困難，已成為化學工程領域基礎問題之一[1]。

由于實驗測量難度的限制，以往的實驗研究大多集中于邊界層中物質擴散現象，而對于湍流中傳質擴散過程實驗研究則不多，實驗手段也較少，其中較為有效的包括采用粒子成像測速儀（PIV）針對流場和激光誘導熒光（LIF）針對濃度場測量研究[3-5]，有的結合了激光多普勒測速（LDV）技術[6]。流型包括帶軸向葉輪的攪拌槽[5]以及帶有格柵的水平流體通道[6]。

在湍流中傳質擴散過程數值模擬方面，近年來的研究主要集中在基于嚴格Navier-Stokes 方程的湍流擴散模型研究。例如Sun 等[7]通過湍流傳質模型計算了湍流傳質擴散系數，并將其用于板式精餾塔的數值模擬中；Liu 等[8]通過數值模擬，預測了散堆填料精餾塔中的液相濃度分布，預測結果與實驗測量結果相吻合；Dong 等[9]利用3D 體積流體模型（VOF）建立了對流擴散模型，模擬了二氧化碳在氫氧化鈉溶液中的反應吸收過程，計算出了液體傳質系數；Zhang等[10]將湍流傳質的理論模型應用到鼓泡塔反應器中，同時精確地預測鼓泡塔中物質的濃度和速度分布。雖然上述研究建立了湍流條件下傳質的嚴格模擬方法，但是，因為嚴格模擬依賴于計算微元網格上的離散，導致必須以巨大的計算量來換取計算精度，因此，計算量與計算精度之間權衡難以實現，嚴重制約著這一類方法的應用和發展。因而針對湍流傳質擴散過程建立新的更加快速、高效的數值方法具有重要的意義。

如果將湍流作為一種時間序列隨機過程，那么對其進行數值描述面臨超高的維度。為此研究者提出了湍流流場的降維重構方法，即使用有限的特征來對復雜系統進行近似。這一種降維方法的關鍵是從已知的流場數據中尋找和提取復雜物理過程中主要的特征，或稱模態[11-12]。本征正交分解（proper orthogonal decomposition,POD）是一種常用的降維和特征提取方法，能夠有效獲取湍流過程中的一系列相干結構以及每個結構重要性參數，即能量[13]。Zhang等[14]測量了平板分離再附流動的速度場，并對實驗數據進行了POD 分析，通過低階特征模態重構了流場結構。Lacassagne 等[15]利用粒子成像測速（PIV）技術，測量了剪切變稀聚合物溶液中振蕩格柵湍流的流場，并應用POD對流場進行了周期性分析。Taira 等[16]以數值模擬的數據為基礎，應用POD 方法提取了圓柱繞流流場的空間模態，并分析了圓柱繞流的周期性。然而已有研究都是對速度場的特性進行POD 分析，迄今尚未有關于湍流中的濃度分布進行POD 分析研究的報道。POD 是一種基于數據的復雜過程數值重構方法，建立基于POD 的湍流傳質擴散的數值模擬方法具有重要的探索價值。

本文以熒光素鈉水溶液在格柵通道中的湍流擴散過程為研究主體，采用PIV 和LIF 聯用方法，同時測量了湍流通道的速度和濃度瞬時分布，并對實驗結果進行定量分析，應用POD 方法對中心位置處測得的瞬時濃度分布進行了模態分解及周期性分析。最后，利用POD 分解得到的低階模態，進行了湍流擴散濃度分布的數值重構，為快速預測湍流擴散中的濃度分布提供了基礎。

1 實驗原理和方法

1.1 測量原理

PIV 技術原理是通過CCD 相機，借助于具有良好跟隨性示蹤粒子的散射光，記錄流體某一區域中多個粒子的位置，并計算兩次曝光時間間隔內每個粒子的位移矢量，進而得到該區域速度矢量場[17-18]。由于湍流狀態下流體的速度較大，故本實驗PIV 相機采用雙幀雙曝光模式。

LIF 技術原理是用特定波長的激光照射溶液中熒光劑來激發熒光，通過CCD 相機記錄熒光強度分布，并將光強轉換為濃度進而獲得溶液中熒光劑的濃度分布[19]。當熒光劑濃度較低時，可以忽略激光光束和熒光信號在吸收介質中的衰減，熒光強度和熒光劑濃度之間有如下關系[2]：

式中，If為被檢測的熒光強度；Kopt為光學常數；I0為入射激光強度；φ為熒光量子產率；ε為熒光劑摩爾吸收系數；V為采集體積；Cf為熒光劑的濃度。當熒光劑、裝置位置、激光強度等實驗條件固定時，除熒光劑濃度變量外，等式右邊其他系數都可視為常數，則熒光強度與熒光劑濃度呈線性關系。

1.2 實驗試劑

實驗采用熒光素鈉（C20H10Na2O5）（純度≥98%）作為熒光劑，其極易溶于水且熒光較強的特性適用于LIF 濃度測量；示蹤粒子是空心玻璃球（直徑8～12 μm，密度900 kg·m-3），其與水的密度十分接近，在水中的跟隨性良好，粒子的運動就代表了湍流的流動，因此示蹤粒子對熒光劑擴散的影響可以忽略；水為經凈化處理的自來水。

1.3 實驗裝置

如圖1(a)所示，本文所用實驗裝置由液相湍流擴散發生系統和PIV/LIF 數據采集系統組成。液相湍流擴散發生系統是由儲水槽、潛水泵、流體通道、轉子流量計以及熒光素鈉溶液注入系統組成。其中流體通道為有機玻璃制成的透明方形截面通道，內部尺寸為1000 mm×60 mm×60 mm（長×高×寬），如圖1(b)所示。熒光素鈉溶液通過隔膜泵、轉子流量計和圓形噴嘴組成的注射裝置在通道截面中心位置注入。為了讓熒光劑的擴散在較為均勻（充分發展）的湍流中進行，熒光劑注入的軸向位置選在距通道進口400 mm處，同時在其上游70 mm處設有金屬絲網格柵，絲徑為1 mm，柵格間距為6 mm。

圖1 實驗裝置Fig.1 Schematic diagram of experimental apparatus

聯合PIV/LIF 數據測量系統包括一個Nd∶YAG激光器和兩個CCD 相機，其中Nd∶YAG 激光器輸出波長為532 nm，最大脈沖能量為200 mJ 柱狀激光，通過片狀光學器件將柱狀激光重塑為最小厚度約為1 mm 的楔形片狀激光光束垂直穿過擴散通道。兩個CCD相機分別用于PIV和LIF系統，用來記錄流體中激光照射的片狀區域中示蹤粒子的流動位移和熒光素鈉的熒光強度，便可獲得照射的片狀區域二維信號分布。CCD 相機的分辨率均為1376 pixel×1040 pixel。PIV 相機前加設波長為532 nm 的窄帶濾光片來消除背景光的影響，LIF 相機則加設550 nm 以上的高通濾鏡來保留熒光和濾除激發光。為了同時測量相同位置湍流速度和濃度分布，PIV相機與LIF 相機位于同一高度且聚焦相同的視場，由兩個相機方向之間的夾角導致的圖像誤差由軟件加以校正。兩個相機和激光器的操作控制以及圖像處理和校正通過LaVision 公司的DaVis8.2 軟件同步完成，圖像的拍攝頻率為5 Hz。

1.4 實驗方法

實驗開始前，需要對LIF 熒光強度與實際熒光劑濃度關系[即式（1）中的常數]進行標定[20]。首先配制不同濃度熒光素鈉溶液，再利用LIF 相機分別拍攝50 張連續圖像，去除背景噪聲后取平均值，即得到不同濃度熒光素鈉溶液對應的熒光強度值。得到的標定曲線如圖2 所示。可以看出，熒光強度和熒光劑濃度間基本呈線性關系。

圖2 熒光強度-熒光素鈉濃度標定曲線Fig.2 Calibration curve between fluorescence intensity and C20H10Na2O5concentration

在實驗過程中，先向擴散通道中注入水，控制水的流速為0.2 m·s-1，對應的Reynolds數約為12000，流動狀態為湍流。待流速穩定后，將300 mg·L-1的熒光素鈉溶液注入流體通道，同時保證熒光素鈉溶液的流速也為0.2 m·s-1。通過電腦控制激光器、LIF 相機和PIV 相機進行同時拍攝，記錄100 s時間內熒光素鈉溶液在水中擴散過程。實驗過程中房間應保持黑暗以避免測量噪聲。實驗在（25.3±0.3）℃下進行，擴散裝置、激光頭、相機均固定于防震平臺上。

在PIV圖像處理過程中，判詢域尺寸為48 pixel×48 pixel，利用時間序列互相關算法計算各個判詢域中示蹤粒子的位移，從而計算出速度矢量即為流體的瞬時速度。LIF圖像的處理如圖3所示，首先去除背景噪聲，再使用DaVis8.2 軟件中位移校正功能將處理后的LIF圖像調整到與PIV圖像相同的坐標，利用圖2的標定曲線將LIF圖像熒光強度轉換為熒光素鈉的濃度。實驗中通過改變激光片狀光束與流體通道的相對位置，可以獲得流體中不同位置的二維信號。

圖3 LIF圖像處理流程圖Fig.3 Flow diagram for LIF image processing

2 實驗結果與討論

2.1 湍流擴散濃度和速度分布

應用上述實驗和數據處理方法，可以得到不同時刻t下，中心切面（z=0）上的熒光素鈉濃度分布，結果如圖4 所示。縱坐標y表示觀測區域的高度，y=0處為熒光素鈉噴嘴的位置，橫坐標x表示觀測區域的長度，坐標原點為熒光劑入口。由圖4 熒光劑濃度分布可以看出流體處于較強程度的湍流狀態。

圖4 熒光素鈉在湍流水中擴散過程的不同時刻熒光素鈉濃度分布（中心位置z=0）Fig.4 Concentration distribution during diffusion process at different times(central position z=0)

對100 s 內拍攝的500 張中心切面的LIF/PIV 圖像，分別求時間平均后可以得到平均濃度和平均速度分布，如圖5所示。從圖5(a)可以看出，受軸向流體流動的影響，軸向擴散程度比徑向高，因此熒光素鈉的時均濃度分布呈錐形。從圖5(b)可以看到，平均后的速度矢量基本沿水平方向，并且由于湍流邊界層在壁面上的發展，沿著流體流動方向，壁面附近平均速度逐漸減小，中心區域的平均速度逐漸增加。

2.2 不同切面位置濃度分布

為了進一步探究湍流擴散的特性，本文還測量了不同切面位置處的濃度分布，其時間平均濃度分布如圖6所示。不同切面位置的濃度分布是通過保持激光和CCD 相機的位置不變，然后前后平移流體通道裝置，因介質厚度改變從而微調相機焦距來獲得的，z為偏離中心切面的距離。為了使不同切面的熒光素鈉濃度顯示得更清楚，縮小了平均濃度的標尺，如圖6(a)所示。圖6(b)、(c)是距離中心切面8 mm位置的平均濃度分布，可以看出兩圖的濃度分布形狀基本一致，均呈現錐形分布，并且沿軸向方向0～20 mm 處的熒光素鈉濃度較低，20～100 mm 間濃度較高，而后濃度又減小，說明此擴散過程呈現三維的錐形分布，即入口附近湍流擴散程度較小，隨著軸向距離增加，擴散程度逐漸增大。圖6(d)、(e)為距離中心切面16 mm切面位置的平均濃度分布，與8 mm的濃度分布相比，其靠近入口處的低濃度區域進一步增大，后續錐形區域內的濃度也更低。

圖6 不同切面位置平均濃度分布Fig.6 Mean concentration distribution of different section positions

2.3 濃度分布定量分析

圖7 為沿z方向不同切面的熒光素鈉時均濃度分布，其中橫坐標是y方向距離，縱坐標分別是熒光素鈉時均濃度或歸一化的時均濃度。每一切面都畫出了沿流動方向x=0～50 mm 處的時均濃度變化曲線。圖7(a)為z=0（中心位置）處不同軸向距離上的時均濃度曲線。y=0即熒光劑入口位置，熒光素鈉濃度最高，隨著軸向距離的增加，濃度逐漸減小，由于受湍流擴散的影響，濃度曲線均呈高斯分布。沿流動方向較遠的x=30、40、50 mm處的濃度曲線變化較平緩，為了清楚觀察其濃度分布特征，對圖7(a)的時均濃度進行了歸一化處理，結果如圖7(b)所示。圖7(b)中，Cc表示熒光劑入口中心線(y=0)處的濃度值，歸一化后的濃度同樣符合高斯分布。圖7(c)、(d)分別為z=8 mm 和z=-8 mm 切面上不同軸向時均濃度曲線，可以看出兩圖曲線變化基本一致，說明湍流擴散的時均濃度分布呈軸對稱結構。此外，不同軸向濃度變化與z=0相反，隨著軸向距離增加，濃度逐漸增大，這是因為湍流擴散為錐形即擴散程度在逐漸增大。圖7(e)、(f)分別為z=16 mm 和z=-16 mm 對應的平均濃度曲線，其變化規律與圖7(c)、(d)相似。

圖7 不同切面平均濃度分布軸向方向距離入口0～50 mm處平均濃度變化Fig.7 The mean concentration distribution of different sections at 0 to 50 mm away from the outlet in the axial direction

3 濃度分布POD模態分析

3.1 POD模態分解原理

POD 方法是一種高效的數據降維分析方法，其目標是把復雜的高維物理過程用低維模態進行近似，進而顯著減小用數值方法重現復雜物理過程所需的計算量和存儲數據量[21]。湍流并非是完全不規則的隨機運動，而是存在一定結構特征的復雜流體運動，在這樣的流場中，熒光劑濃度作為一個物理參數，其分布存在一系列可識別的具有周期性的結構，這種有序結構稱為相干結構[22]。而通過POD分解所得到的模態可以識別出濃度分布的相干結構[23]。

POD 方法將瞬時濃度矩陣C(x)分解為時間系數{θi}和空間模態{φi(x)}[24-27]序列，其中i=1,…,N，N為完全重構瞬時濃度分布所需的模態數，即

瞬時濃度數據矩陣由LIF實驗獲得，假設在LIF圖像上有m個離散采樣點，單次采樣就可以同時得到m個離散點的值，共計采集N次。以矩陣的形式表示為：

POD 算法可以求出一組模態[28-29]，通過這組模態可以得到瞬時濃度分布的最佳近似值，同時每個模態滿足與所有其他模態正交的條件。為建立等效的特征值問題，由C建立協方差矩陣R

R為N×N階矩陣，其維數遠低于C。于是空間模態可表示為：

式中，λ是模態φ對應的特征值，表示相應模態的能量大小，并由大到小進行排列[30]。

時間系數矩陣θ由式（6）計算得到：

3.2 濃度分布的POD模態分解

從圖4所示的瞬時濃度分布，可以看到在遠離熒光劑入口區域的濃度很低且接近于0，而對這些區域進行瞬態POD 模態分析，將導致大量的數據和計算資源被用來處理這些不明顯區域。而且，低濃度區域的濃度受測量噪聲的影響更大，而這些噪聲是無法有效進行POD 模態分析的。因此，本文僅對熒光劑入口區域附近進行了POD模態分析。截取的瞬時濃度分布如圖8所示，x=0處為熒光劑入口位置。

圖8 截取后的瞬時濃度分布（t=10 s）Fig.8 Instantaneous concentration field after interception(t=10 s)

應用式（6）得到湍流擴散的瞬時濃度分布的前四階模態，如圖9 所示。從圖9(a)可以看出，1 階模態中，正負時間系數呈現上下對稱的分布結構；圖9(c)、(d)所示的3 和4 階模態大致在0～10 mm 范圍內呈現前后交替的結構；圖9(b)為2 階模態，處于過渡階段，既有前后交替又有上下反對稱結構。從模態分解的結果可以定性看出，在x方向前10 mm 的區間內，湍流擴散的濃度分布具有較明顯的規則結構特征，10 mm以后則較弱。

圖9 POD前4階模態θ值的分布Fig.9 Distribution of θ of the first 4 POD modes

3.3 濃度分布模態結構周期性分析

基于圖9 所示的結構分布特征，可對實驗區域進行分區，現將距離熒光劑入口較近的x=0～10 mm的區域稱為A 區，較遠的即x＞10 mm 區域稱為B 區。前4階模態沿軸向的平均時間系數分布曲線如圖10所示，可以看出，前4 階模態的平均時間系數在A 區規律性更強，均呈振幅不同的類正弦曲線分布，而各階模態在B 區的曲線則更不規則。在各階模態中，類似正弦的時間系數曲線說明此區域的瞬時濃度分布周期性更明顯[17]。因此，在A 區，湍流擴散導致的熒光劑濃度分布具有明顯的周期性，隨流動距離增加，湍流擴散的周期性逐漸減弱、消失。

圖10 平均時間系數Fig.10 Mean temporal coefficient

3.4 濃度分布POD模態重構

從模態分解和分析結果可以看出，湍流擴散在模態空間中具有一定的規律。因此，基于以上分解結果，可以對湍流擴散的濃度分布進行重構。在式（5）中，特征值可以表示相應模態的能量大小，模態的能量越大，表示該階模態對湍流濃度分布的貢獻度和影響越大。并且POD分解得到的模態是按其能量貢獻大小進行降序排列的[31],因此低階模態占有更高的能量比例。各階模態的能量占總能量的比例如式（7）所示，式中λi為第i階模態對應的特征值[32]。前p階模態的能量占總能量的比例如式（8）所示。

圖11(a)為A 和B 區前500 階POD 模態能量圖，可以看出在A 和B 區內，低階模態均含有較高的能量，但A 區的低階模態能量占比高于B 區，說明A 區的能量在低階模態更集中。隨著模態階數的增加，能量占比急劇下降，并逐漸接近于0。圖11(b)為前500 階模態的能量占比的累積值，可以發現A 區前100階模態的能量之和占前500階總能量的85%，而B區前100階模態的能量僅占68%，表明周期性較強的A區應用POD模態重構會更為有效。

圖11 POD模態能量（A和B區）Fig.11 POD modal energy(A and B regions)

根據式（2），由前r階模態所重構的濃度分布為[33]：

圖11(b)表明，隨著r的增加，重構的精度隨之增加，但會顯著減緩。

圖12 是分別選擇前10、20、50 和100 階模態重構得到的瞬時濃度分布。與圖8 比較可以看出，選用的階數越高，重構得到的濃度分布越接近于實驗濃度分布。但隨著階數的增加，計算量和所需要的數據量也隨之增加。圖12 利用前100 階模態重構的濃度分布基本反映了原始分布的主要特征。

圖12 POD模態重構Fig.12 POD modal reconstruction

圖13 為在x方向上不同位置重構的濃度分布與實驗值的比較結果。可以看出，當x≤10 mm，50和100 階重構結果均能較好地捕捉實驗濃度分布特征；x=15 mm 處僅前100 階重構較為接近實驗結果；而x=20、25 mm 的重構結果均與原始相差較大，表明A 區低階和高階重構效果均較好，超出該區域則需要較高階模態重構，以致不能準確重構。

圖13 不同軸向位置的模態重構和原始濃度分布Fig.13 Modal reconstruction and original concentration distribution of different coaxial positions

圖14 給出了4 種不同階數重構結果與實驗值之間的相對誤差分布。可以看出，不同階數的重構結果中，誤差較大的區域均集中于B 區右下方，且隨著階數的升高，高誤差區域逐漸減小，這可能是因為湍流在后端逐漸發展，速度邊界層逐漸建立，內部流動逐漸接近橢圓側型流動，此時的流動狀態與前端大不相同，具有某種傾向性，壓抑了隨機性，加之邊壁效應，使得誤差增大。而A 區的相對誤差均在1% 以內，重構結果十分接近實驗結果。

圖14 模態重構相對誤差Fig.14 Relative error of modal reconstruction

對上述相對誤差圖進行分區處理，分別對A和B 區的相對誤差求平均，得到如圖15 所示曲線。可以發現，兩個區域的平均相對誤差均隨模態階數的增加而減小，A 區的平均相對誤差極低，前10階模態重構的平均相對誤差為0.74%，前100 階模態重構的平均相對誤差僅為0.45%；而B 區的平均相對誤差相對較高，前10 階模態重構的平均相對誤差為11%，前100 階模態重構的平均相對誤差為3.8%。誤差的分析結果進一步說明了具有周期性的A 區域POD 模態重構效果更好，重構誤差集中在B 區。

圖15 模態重構平均相對誤差（A和B區）Fig.15 Mean relative error of modal reconstruction(A and B regions)

4 結論

（1）本文采用PIV 與LIF 聯用技術實現了湍流擴散過程速度和濃度分布的同時測量，發現此湍流擴散過程的時間平均濃度分布呈三維錐形，且不同切面位置的徑向平均濃度均呈高斯分布。

（2）建立了利用POD 模態分解方法對瞬時濃度分布進行數值重構的方法，將濃度分布進行模態分解，得到有一定規律的時間系數分布。

（3）在接近熒光劑入口區域內，低階模態的平均時間系數呈現振幅不同的類正弦曲線分布，具有周期性，而遠離熒光劑入口區域的模態沒有此特征；接近熒光劑入口區域的低階模態時間系數在模態空間內的分布呈現一定的相干結構和規律，而遠離熒光劑入口區域則沒有明顯規律。

（4）接近熒光劑入口區域和遠離該區域的各階模態能量占比變化具有相似趨勢，即低階模態均具有較高能量，隨著模態數的增加，能量占比急劇降低，逐漸接近于0。但接近熒光劑入口區域的能量更集中在低階模態，前100 階模態能量占85%，而遠離熒光劑入口區域前100階模態能量僅占68%。

（5）對瞬時濃度分布進行了POD 模態重構及誤差分析，越高階數重構所得到的濃度分布越接近于原始濃度分布，且包含的主要特征越多；近熒光劑入口區域和遠離該區域的平均相對誤差均隨重構模態數的增加而減小，且重構誤差主要集中在遠離熒光劑入口區域，結果表明周期性較強的近熒光劑入口區域更適合應用POD 模態進行濃度分布的重構。考慮到計算量應采用較少的模態數，在較為規則的湍流區域，使用POD 較低階模態就可以預測，而用普通的直接數值模擬難以奏效。