分形流場傳質(zhì)輸運特性的數(shù)值模擬研究

石 均，陶善聰，周 毅

(1.南京理工大學(xué)能源與動力工程學(xué)院，江蘇 南京 210094；2.中國航天空氣動力技術(shù)研究院，北京 100074)

湍流現(xiàn)象廣泛存在于自然界和工業(yè)界中，如森林火災(zāi)的大規(guī)模傳播[1，2]、城市熱島效應(yīng)[3]、發(fā)動機(jī)燃燒室的湍流燃燒[4]、機(jī)翼尾流[5]等.根據(jù)已有研究表明，湍流可以提高流動系統(tǒng)的輸運能力和擴(kuò)散混合強(qiáng)度[6].近些年，眾多學(xué)者采用不同的方式來研究各種湍流問題.其中較為典型的是利用仿生學(xué)思想[7]，提出將分形理論[8]與湍流領(lǐng)域相結(jié)合，從而來探索更為復(fù)雜的湍流現(xiàn)象，以期為實現(xiàn)相關(guān)的工業(yè)應(yīng)用提供理論指導(dǎo).

分形網(wǎng)格湍流FGT(Fractal-Generated Turbulence)是流體流過具有分形特征的格柵時引發(fā)的不規(guī)則、混亂的流動現(xiàn)象，如上文提到森林火災(zāi)的擴(kuò)散呈現(xiàn)出分形湍流的特點[1，2].與之對應(yīng)的則被稱為規(guī)則網(wǎng)格湍流RGT(Regular-Generated Turbulence)，即流體流經(jīng)具有規(guī)則結(jié)構(gòu)的格柵后產(chǎn)生的湍流現(xiàn)象.此前，經(jīng)研究發(fā)現(xiàn)分形流場具有比規(guī)則流場更強(qiáng)的湍流擴(kuò)散混合能力[9]，為此人們對這種分形湍流過程的研究熱度與日俱增.Queiros-Conde等[10]早期便是利用在風(fēng)洞中加入分形格柵來開展分形流動研究，結(jié)果發(fā)現(xiàn)分形流動的結(jié)構(gòu)函數(shù)尺度指數(shù)變化趨勢與分形維度大小相反.Staicu等[11]利用分形流動實驗探索了慣性尺度和耗散尺度范圍內(nèi)的流場特點，證明分形湍流的自相似性行為受到分形結(jié)構(gòu)維度的影響.Hurst和Vassilicos[12]開展分形湍流實驗，對比了十字形、工字形和正方形三種網(wǎng)格形成的湍流場特性，結(jié)果發(fā)現(xiàn)流場下游的均勻各向同性取決于分形結(jié)構(gòu)的幾何維度或網(wǎng)格的尺寸.Seoud和Vassilicos[13]對上述幾種網(wǎng)格湍流進(jìn)行了更深入的研究后發(fā)現(xiàn)，在遠(yuǎn)離網(wǎng)格的流場下游衰減區(qū)域，泰勒微尺度與積分尺度的變化保持相對獨立，并與入口初始速度無關(guān).Laizet和Vassilicos[14]采用直接數(shù)值模擬DNS(Direct Numerical Simulation)手段對FGT和RGT進(jìn)行仿真研究，發(fā)現(xiàn)分形流場具有更高的渦度和湍流強(qiáng)度.Laizet等[15，16]利用高保真數(shù)值方法Incompact 3D對FGT進(jìn)行一系列研究后指出，在分形流場遠(yuǎn)下游區(qū)域的小尺度尾流會被隨機(jī)卷入大尺度尾流中，并且不停地重復(fù)這種行為.因此他們認(rèn)為是特殊的分形結(jié)構(gòu)導(dǎo)致了分形流場的標(biāo)量輸運和擴(kuò)散能力更強(qiáng)，同時據(jù)此提出了一種空間尺度展開SSU(Space-Scale Unfolding)機(jī)制.此后，Laizet和Vassilicos[17]開展了二維分形湍流的數(shù)值模擬，證實了由于SSU機(jī)制導(dǎo)致分形流場的標(biāo)量橫向擴(kuò)散能力大于規(guī)則流場.另外，特別是近年來，大量學(xué)者積極探索分形湍流理論的其他應(yīng)用可能性.例如，Cafiero等[18，19]基于分形網(wǎng)格增強(qiáng)射流湍流強(qiáng)度的原理進(jìn)行實驗研究，提出了一種圓形沖擊射流被動強(qiáng)化傳熱的新方法.Paul等[20]對單一方形結(jié)構(gòu)湍流場開展DNS研究，通過對Nusselt的統(tǒng)計分析發(fā)現(xiàn)，與前沿駐點相比，流場最大熱傳遞點附近發(fā)生了極端事件.并且他們建議可以將相應(yīng)的研究結(jié)果應(yīng)用于換熱器傳熱優(yōu)化等實際問題.Hoi等[21，22]將分形結(jié)構(gòu)應(yīng)用于翅片散熱器中，在換熱器領(lǐng)域取得了大量研究成果.Gehlert等[23]將分形端板引入到NACA0012半翼展中進(jìn)行飛行實驗，結(jié)果發(fā)現(xiàn)分形端板對翼尖渦結(jié)構(gòu)和尾流場特性產(chǎn)生了較大影響，這意味著可將分形結(jié)構(gòu)應(yīng)用于飛行器部件設(shè)計中.

此外，分形結(jié)構(gòu)提高湍流擴(kuò)散混合能力的特點在諸多工業(yè)設(shè)備中具有廣泛應(yīng)用潛力，如采用分形刀刃的攪拌器[24]及燃燒設(shè)備中設(shè)置分形結(jié)構(gòu)[25]等.然而，經(jīng)過文獻(xiàn)分析發(fā)現(xiàn)，大量學(xué)者關(guān)于分形流場的研究大都集中探索湍流統(tǒng)計特性的機(jī)理方面，尚未有學(xué)者對分形流場的傳質(zhì)輸運特性進(jìn)行深入分析.因此，本研究在課題組前期工作[26]的基礎(chǔ)上，擬針對分形流場的輸運機(jī)理開展研究.本文基于偽譜Fourier-Galerkin方法，將分形湍流按時間演化方式進(jìn)行直接數(shù)值模擬.當(dāng)分形流場演化至一定狀態(tài)后，采用拉格朗日追蹤方法，研究流場中流體微元的運動規(guī)律.同時模擬了具有相同阻塞率的RGT，以作對比分析.

1 直接數(shù)值模擬的方法簡介

1.1 數(shù)值模擬的控制方程與數(shù)值方法

本研究的DNS仿真采用Hussaini等[27]提出的偽譜方法，流場在物理空間的控制方程包括連續(xù)性方程和Navier-Stokes方程(簡稱NS方程)，如公式(1)和公式(2)，其中空間網(wǎng)格節(jié)點進(jìn)行均勻離散.

?·u=0，

(1)

(2)

公式中：u為瞬時速度矢量；ν為運動粘度；P=p+u·u/2為無量綱化后的壓力值.

在偽譜伽遼金方法中，需要將變量進(jìn)行離散傅里葉變換至譜空間求解，因此控制方程最終在譜空間的為公式(3)和公式(4)，具體的求解過程見參考文獻(xiàn)[26].

(3)

(4)

本研究數(shù)值方法的時間遞進(jìn)采用四階龍格-庫塔格式.根據(jù)顯示格式要求庫朗數(shù)CFL(Courant Friedrichs Lewy)小于1，設(shè)置時間步長為Δt=0.005M/U0(M為規(guī)則相鄰網(wǎng)格條之間的間距，同時也等于分形網(wǎng)格第三階正方形的長度尺寸；U0為流場流向的初始平均速度).此外，NS方程的對流項、粘性項和壓力項均采用譜方法求解.

1.2 計算域的參數(shù)設(shè)置

圖1 分形網(wǎng)格與規(guī)則網(wǎng)格Y-Z平面示意圖

此外，下面將給出數(shù)值模擬的部分參數(shù)，具體如表1所示.

表1 計算工況參數(shù)

兩類流場網(wǎng)格在Y-Z平面上用Nf表示幾何維度，規(guī)則網(wǎng)格的Nf=1，而分形網(wǎng)格的Nf=3(分別對應(yīng)i=0，1，2).網(wǎng)格阻塞率為σ=Tbar/(LYLZ)，Tbar表示Y-Z面網(wǎng)格條面積，LY和LZ分別為流場沿法向和展向的長度.本研究流場計算域的長度Lx=LY=LZ=8M.雷諾數(shù)設(shè)置為1 600，ReM=(U0M)/ν基于初始平均速度U0和網(wǎng)格尺寸大小M.將最大正方形網(wǎng)格的橫向厚度與最小正方形橫向厚度的比值定義為分形網(wǎng)格橫向厚度比Dr，即Dr=d0/d2=9.5.

圖2 加入附加擾動后的流場初始狀態(tài)

此外，沿著流向X、法向Y和展向Z的網(wǎng)格節(jié)點數(shù)設(shè)置為Nx×NY×NZ=800×400×400，三個方向的節(jié)點為均勻布置，同時在X、Y和Z方向均采用周期性邊界條件.其余具體的流場設(shè)置細(xì)節(jié)可參考文獻(xiàn)[26]，同時該文章中已經(jīng)給出了FGT和RGT時間演化的初始速度條件.由于網(wǎng)格湍流的時間演化過程中需要向初始流場加入附加擾動，才能加速網(wǎng)格條產(chǎn)生的尾流從層流演變?yōu)橥牧?本文FGT和RGT初始流場中添加的附加擾動場的脈動速度值大約為0.5%U0，特征長度尺度約為0.2M.將此附加擾動速度加入流場后，形成了兩類網(wǎng)格湍流的初始狀態(tài)，如圖2所示.此外，關(guān)于數(shù)值方法準(zhǔn)確性驗證等方面的內(nèi)容詳情亦可參考文獻(xiàn)[26]，在此不再贅述.

2 粒子追蹤期間的流場特性簡介

流場演化至一定狀態(tài)后，選取統(tǒng)計時間段為0≤Ttrack≤100Tr，主要研究流體微元的運動規(guī)律，下標(biāo)“track”代表追蹤粒子的含義.為了便于表述，后文統(tǒng)一將統(tǒng)計示蹤粒子運動的初始時刻定義為Ttrack/Tr=0，結(jié)束統(tǒng)計的時刻定義為Ttrack/Tr=100，Tr=M/U0為無量綱化時間.需要明確的是，其中統(tǒng)計示蹤粒子運動的初始時刻Ttrack/Tr=0并不等同于圖2所示的流場初始演化時刻.因此，接下來首先對0≤Ttrack≤100Tr范圍內(nèi)的流場進(jìn)行初步的介紹，旨在更清晰認(rèn)識此時間段內(nèi)流場的具體特征.

圖3 不同統(tǒng)計時刻下FGT的流向平均速度云圖

2.1 流向平均速度場

圖3和圖4分別為FGT和RGT在不同統(tǒng)計時刻下的流向平均速度〈U〉x云圖，選取了Ttrack/Tr=0和100時刻的結(jié)果進(jìn)行展示，符號“〈?〉x”表示對X方向進(jìn)行平均.對于圖3(a)的FGT而言，當(dāng)Ttrack/Tr=0時，盡管分形流場Y-Z平面上整體的流向平均速度滿足〈U〉x/U0≈0.36，但明顯存在大尺度結(jié)構(gòu)的痕跡，這說明此時的流場依然受到第一階最大分形網(wǎng)格的影響.隨著時間的演化，在圖3(b)的Ttrack/Tr=100時，該大尺度結(jié)構(gòu)的痕跡逐漸弱化，但其并非位于流場的中心區(qū)域，而是運動到了靠近流場的邊界處.針對圖4的規(guī)則流場演化結(jié)果，可以看出Ttrack/Tr=0至Ttrack/Tr=100，規(guī)則流場中幾乎不存在大尺度痕跡，而且整個Y-Z平面的流向平均速度滿足x/U0≈0.36，達(dá)到了較為均勻的狀態(tài).所以對比圖3和圖4可知，在統(tǒng)計時間(0≤Ttrack≤100Tr)內(nèi)，分形條件依然影響著流場的發(fā)展?fàn)顟B(tài).

圖4 不同統(tǒng)計時刻下RGT流向平均速度云圖

圖5 流場的統(tǒng)計特性隨時間演化結(jié)果

圖5(a)為流場中心線上的流向平均速度Uc隨時間演化結(jié)果，統(tǒng)計時間段內(nèi)的FGT和RGT均滿足Uc/U0≈0.36，因此流場沿流向發(fā)展基本趨于統(tǒng)計均勻，其中下標(biāo)“c”代表流場中心線.根據(jù)課題組前期文章[26]的研究可知，初始條件對流場的影響消失后，整個流場都將演化至近似均勻各向同性狀態(tài).圖5(b)為統(tǒng)計時間(0≤Ttrack≤100Tr)內(nèi)流場展向與法向脈動速度比值，可以發(fā)現(xiàn)FGT和RGT兩者的〈w′〉1/2/〈v′〉1/2均在0.85～1.25之間波動，符號“〈?〉”表示對X、Y和Z三個方向進(jìn)行空間平均.其中，F(xiàn)GT的〈w′〉1/2/〈v′〉1/2比值從1.1開始隨時間減小，而RGT的〈w′〉1/2/〈v′〉1/2在1.2附近徘徊.Laizet和Vassilicos[14]進(jìn)行的FGT和RGT空間演化研究中也觀察到，具有相同阻塞率σ時，分形流場和規(guī)則流場下游位置的〈w′〉1/2/〈v′〉1/2在0.8～1.2范圍內(nèi)波動.同時結(jié)合Hurst和Vassilicos[12]的研究結(jié)論可知，該脈動速度比值越接近1，說明此時流場的大尺度各向同性效果越好，并且反映出了流場在Y-Z方向上的對稱性.

2.2 Q值瞬時等值面

為了進(jìn)一步探討兩類流場在統(tǒng)計時間(0≤Ttrack≤100Tr)內(nèi)擬序結(jié)構(gòu)發(fā)展的異同，采用速度梯度張量第二不變量Q值來進(jìn)行識別，具體結(jié)果如6和圖7所示.其中速度梯度的張量Aij=?ui/?xj包括兩部分，瞬時應(yīng)變率張量sij=(?ui/?xj+?uj/?xi)/2和瞬時旋轉(zhuǎn)率張量ωij=(?ui/?xj-?uj/?xi)/2，根據(jù)定義，速度梯度張量第二不變量Q=Qω+Qs=(ωijωij/2-sijsij/2).FGT選取的Q值瞬時等值面為Q/(1/Tr2)=4×10-4.根據(jù)之前的文章[26]研究可知，RGT的衰減速度要快于FGT，其渦度值也更低.為了在相同的統(tǒng)計時刻下清晰識別RGT的擬序結(jié)構(gòu)，RGT選取的Q值瞬時等值面為-1×10-5.

圖6 FGT不同時刻下的Q值瞬時等值面

在圖6和圖7中，左圖是流場擬序結(jié)構(gòu)三維效果，右圖為流向位置X=4M的平面上擬序結(jié)構(gòu)二維效果.在圖6的分形流場中，存在尺度大小不一的擬序結(jié)構(gòu)，但明顯以大尺度結(jié)構(gòu)為主.而且大尺度結(jié)構(gòu)的位置主要靠近邊界區(qū)域，外形類似于經(jīng)典的“蟲子結(jié)構(gòu)”[28，29].同時圖6(a)和圖6(b)的右圖X=4M二維平面上的擬序結(jié)構(gòu)也可以看出，此過程中大尺度結(jié)構(gòu)普遍在流場邊界區(qū)域運動，這與圖3流向平均速度演化觀察到的現(xiàn)象一致.

圖7 RGT不同時刻下的Q值瞬時等值面

根據(jù)圖7(a)和圖7(b)的規(guī)則流場三維擬序結(jié)構(gòu)發(fā)現(xiàn)，規(guī)則流場在此演化過程中不僅渦度值更小，整體的擬序結(jié)構(gòu)尺度大小并未表現(xiàn)出強(qiáng)烈的反差，反而是擬序結(jié)構(gòu)的尺度相比于分形流場分布顯得更加大小均勻.當(dāng)流向位置X=4M時，該二維Y-Z平面上的大部分?jǐn)M序結(jié)構(gòu)大小基本相同，結(jié)合圖4的Y-Z平面上流向平均速度可視化結(jié)果，說明此時的規(guī)則流場演化至一種更加均勻的狀態(tài).

3 計算結(jié)果與分析

本研究將采用Yeung和Pope[30]提出的13點插值方法來獲取粒子統(tǒng)計單元的拉格朗日數(shù)據(jù)，該方法也被稱為TS13插值方案，并且其有效性也得到了Yeung和Pope[31]的證明.數(shù)值模擬時間步長為Δt/Tr=0.005，在課題組前期文章[26]的流場基礎(chǔ)上，選取20 000個時間步長Δt/Tr，即統(tǒng)計時間段為0≤Ttrack≤100Tr.然后對流場中的流體粒子進(jìn)行隨機(jī)標(biāo)記，追蹤粒子的運動軌跡，從而探索粒子的運動規(guī)律.需要再次強(qiáng)調(diào)的是，此時的統(tǒng)計初始時刻Ttrack/Tr=0并非圖2所示的流場演化初始時刻T/Tr=0，而是流場演化至一定狀態(tài)后，選取其中的某一個時間段進(jìn)行粒子追蹤統(tǒng)計.

3.1 示蹤粒子取樣間隔的無關(guān)性驗證

在Yeung和Pope[31]對TS13插值方法的應(yīng)用中解釋道，需要對示蹤粒子在每一個時間步長下的具體位置和速度進(jìn)行插值計算.由于采用DNS仿真，時間步長為Δt=0.005Tr，本研究在0≤Ttrack≤100Tr輸出的數(shù)據(jù)大約將占用82 TB字節(jié)內(nèi)存，這顯然過于龐大，并且會嚴(yán)重延緩計算速度.根據(jù)2.1和2.2小節(jié)的結(jié)果可知，F(xiàn)GT和RGT兩者的流場在0≤Ttrack≤100Tr內(nèi)發(fā)展至近似均勻狀態(tài)，因此將示蹤粒子的取樣間隔設(shè)置為Δt1=0.05Tr，即每10個時間步長進(jìn)行一次統(tǒng)計.為了驗證取樣間隔對跟蹤粒子計算結(jié)果產(chǎn)生的影響，分別對取樣間隔為0.005Tr和0.05Tr進(jìn)行初步試算，總計試算100個時間步長.在對示蹤粒子軌跡的描述中，主要統(tǒng)計粒子沿流向X、法向Y和展向Z運動的平均位移〈X+〉i/M、〈Y+〉i/M和〈Z+〉i/M以及路程〈SX〉i/M、〈SY〉i/M和〈SZ〉i/M.其中符號“〈?〉i”表示對i個粒子進(jìn)行整體平均.

當(dāng)Ttrack/Tr=0的初始時刻，在分形流場和規(guī)則流場中隨機(jī)標(biāo)記Np=104個流體粒子.圖8和圖9分別是試算周期內(nèi)，F(xiàn)GT和RGT中示蹤粒子沿流向、法向和展向運動的平均位移軌跡結(jié)果.可以看出分形流場中，〈X+〉i/M、〈Y+〉i/M和〈Z+〉i/M均在4附近波動，示蹤粒子在不同取樣間隔時，沿流向運動的最大誤差為0.05%，沿法向和展向運動的最大誤差為0.005%和0.002%.而圖9的規(guī)則流場中，〈X+〉i/M、〈Y+〉i/M和〈Z+〉i/M同樣在4附近波動，沿三個方向運動的最大誤差分別是0.04%、0.004%和0.002%.由此可見，釋放104個示蹤粒子能夠充分采集到流場的拉格朗日信息，此時粒子運動采取不同取樣間隔產(chǎn)生的最大誤差僅為0.05%，這完全可以忽略.并且隨著時間演化，流場將變得更加的各向均勻同性，兩者之間的誤差會逐漸減小.這說明在統(tǒng)計時間0≤Ttrack≤100Tr內(nèi)，示蹤粒子取樣間隔產(chǎn)生的影響可以忽略不計.因此，后文統(tǒng)一將取樣間隔設(shè)置為0.05Tr進(jìn)行計算.

圖8 FGT中不同取樣間隔計算的示蹤粒子運動軌跡

圖9 RGT中不同取樣間隔計算的示蹤粒子運動軌跡

3.2 粒子數(shù)量對擴(kuò)散規(guī)律的影響

由于3.1小節(jié)證明了示蹤粒子取樣間隔對粒子擴(kuò)散無影響，本小節(jié)首先將在流場內(nèi)隨機(jī)標(biāo)記Np=100、1 000、10 000個流體粒子，從Ttrack/Tr=0開始追蹤粒子的運動軌跡，直到Ttrack/Tr=100統(tǒng)計結(jié)束.為了讓讀者對示蹤粒子在流場中的分布情況產(chǎn)生直觀的認(rèn)識，當(dāng)Ttrack/Tr=0，在流場范圍0≤X，Y，Z≤8M內(nèi)隨機(jī)標(biāo)記104個流體粒子，具體效果如圖10所示，其中紅色線條表示粒子初始釋放范圍的輪廓.此外，根據(jù)不同示蹤粒子數(shù)量進(jìn)行平均處理后，可以觀察統(tǒng)計粒子的數(shù)量對流體微元擴(kuò)散規(guī)律的影響.

圖10 Ttrack=0時，104個示蹤粒子在全流場中的分布情況

圖11展示了在分形流場和規(guī)則流場中標(biāo)記不同數(shù)量粒子后，示蹤粒子沿流場流向運動的平均位移軌跡變化結(jié)果.可以看出，無論是分形流場還是規(guī)則流場，示蹤粒子沿流向運動位移〈X+〉i/M始終在4附近成周期性波動.尤其是Np=100時，粒子在0≤Ttrack≤100Tr內(nèi)一共運動了約4.5個周期.由于前文分析可知，該時間段內(nèi)，流場流向平均速度滿足〈U〉x/U0≈0.36，流向計算域長度為LX=8M，經(jīng)計算發(fā)現(xiàn)粒子沿流向應(yīng)該運動約4.5個周期，這與圖11的結(jié)果一致.同時值得注意的是，標(biāo)記的粒子數(shù)量越少，流體粒子更容易受到流場局部大尺度的影響，結(jié)合流場設(shè)置的周期性邊界條件可知，粒子沿著流向運動平均位移將呈現(xiàn)出明顯的周期性波動.然而隨著示蹤粒子數(shù)量增加，可以反映出對流場整體的系綜平均效果.同時結(jié)合系統(tǒng)的質(zhì)量守恒原則，當(dāng)統(tǒng)計全流場中足夠多的示蹤粒子運動時，流體粒子沿流向平均運動位移值可以進(jìn)行理論預(yù)測，即〈X+〉10 000/M=4，粒子沿流向運動平均位移波動幅度越小.這說明隨著統(tǒng)計流體粒子數(shù)量越多，越能夠充分反應(yīng)出流場的真實情況.因此，當(dāng)示蹤粒子數(shù)量Np=104時，已經(jīng)能夠準(zhǔn)確捕捉粒子的運動軌跡，證明粒子數(shù)量此時將不再影響粒子整體的擴(kuò)散效果.

圖11 不同數(shù)量的示蹤粒子沿流向運動位移軌跡

3.3 分形結(jié)構(gòu)對擴(kuò)散規(guī)律的影響

由于課題組前期文章[26]中主要研究了分形流場的統(tǒng)計特性和能量耗散規(guī)律，因此本研究在此基礎(chǔ)上進(jìn)一步探索分形流場的輸運特性.為了觀察分形結(jié)構(gòu)對流體粒子運動規(guī)律的影響，本節(jié)將在分形流場兩個不同的區(qū)域(Run1：中心區(qū)域和Run2：對角區(qū)域)隨機(jī)標(biāo)記流體粒子，從Ttrack/Tr=0開始追蹤粒子的運動軌跡.與此同時，為了進(jìn)行對比，同樣將在規(guī)則流場的對應(yīng)區(qū)域標(biāo)記示蹤粒子，觀察規(guī)則流場中粒子擴(kuò)散規(guī)律與分形流場的差異.Run1的具體詳情為Ttrack/Tr=0時，被動示蹤粒子數(shù)量Np=104，釋放范圍為0≤X≤8M，2M≤Y，Z≤6M.Run1指定的中心區(qū)域以第一階網(wǎng)格產(chǎn)生的大尺度運動主導(dǎo)，且小尺度運動始終被大尺度包圍在中心區(qū)域，該區(qū)域的運動特性最為特殊.另外，Run2的具體詳情為：Ttrack/Tr=0時，被動示蹤粒子數(shù)量Np=104，釋放范圍為0≤X≤8M，4M≤Y，Z≤8M.該對角區(qū)域主要是第二階和第三階網(wǎng)格尺度較為集中，且位置遠(yuǎn)離流場的中心區(qū)域，研究該區(qū)域可探索分形流場中心區(qū)域的大尺度是否會對對角區(qū)域的流體運動產(chǎn)生影響.Run1和Run2指定區(qū)域的網(wǎng)格特點可見圖1所示的網(wǎng)格二維幾何結(jié)構(gòu)，接下來將分別對Run1和Run2的計算結(jié)果進(jìn)行分析.

3.3.1 Run1中心區(qū)域的結(jié)果分析

當(dāng)計算工況為Run1時，統(tǒng)計初始時刻，示蹤粒子在流場中的分布情況如圖12所示.其中的紅色矩形并非真實存在，而是為了便于識別初始的粒子位置所標(biāo)記.圖13為Run1的示蹤粒子在分形流場和規(guī)則流場中運動軌跡隨時間變化效果.在前文已經(jīng)提到，主要統(tǒng)計示蹤粒子沿流場流向、法向和展向運動的平均位移和平均路程.

圖12 Ttrack=0時，Run1的104個示蹤粒子在流場中的分布情況

從圖13(a)的示蹤粒子沿流向運動的平均位移可知，無論是在分形流場或規(guī)則流場中，〈X+〉10 000/M始終在4附近波動，且波動幅度大致相同，說明兩類流場中心區(qū)域沿流向上的湍流特性具有相似性.圖13(b)為示蹤粒子分別沿法/展向運動的平均位移軌跡，對比FGT和RGT發(fā)現(xiàn)，規(guī)則流場中的粒子沿法向或展向運動位移波動幅度較小，〈Y+〉i/M和〈Z+〉i/M幾乎約等于4.而分形流場的粒子沿法向或展向運動時，〈Y+〉i/M和〈Z+〉i/M均呈現(xiàn)出遠(yuǎn)離4的趨勢，說明分形流場中粒子沿法/展向擴(kuò)散強(qiáng)度更大.對于圖13(c)的粒子沿流向運動平均路程而言，F(xiàn)GT和RGT中的示蹤粒子沿流向運動路程與統(tǒng)計時間成線性關(guān)系〈SX〉i≈0.36ΔTtrack，同時Ttrack=100Tr時，兩者的路程均滿足〈SX〉10 000/M≈36.值

圖13 Run1的示蹤粒子運動軌跡

得關(guān)注的是圖13(d)的粒子沿法/展向運動平均路程變化結(jié)果，具體在FGT中，示蹤粒子沿展向運動路程約為法向1.1倍，且當(dāng)Ttrack=100Tr時，〈SY〉10 000和〈SZ〉10 000大約在0.5M～0.6M之間.而在RGT中，示蹤粒子沿展向運動路程約為法向的1.2倍，且當(dāng)Ttrack=100Tr時，〈SY〉10 000和〈SZ〉10 000大約在0.15M～0.2M之間.這說明，兩類流場中的粒子沿展向運動能力都略強(qiáng)于沿法向運動，同時FGT流場的展向和法向各向同性效果更好，這一點從圖5(b)的流場展向與法向脈動速度比值結(jié)果也可以得到證明.此外，F(xiàn)GT中粒子沿法/展向運動路程大約為RGT的3倍，意味著分形流場的示蹤粒子沿法/展向的擴(kuò)散強(qiáng)度比規(guī)則流場更大.

3.3.2 Run2對角區(qū)域的結(jié)果分析

當(dāng)計算工況為Run2時，統(tǒng)計初始時刻，示蹤粒子在流場中的分布情況如圖14所示.圖中的紅色矩形框同樣是為了便于識別初始的粒子位置所標(biāo)記.分形流場的該區(qū)域以第二階和第三階網(wǎng)格尺度尾流運動為主，同時位于邊界處的流體粒子沿法向或展向運動會不停地穿越邊界，結(jié)果可與圖12的結(jié)果進(jìn)行對比.

圖14 Ttrack=0時，Run2的104個示蹤粒子在流場中的分布情況

圖15 Run2的示蹤粒子運動軌跡

Run2工況的計算結(jié)果見圖15.圖15(a)與圖13(a)一樣，F(xiàn)GT和RGT中示蹤粒子沿流向運動的平均位移始終在4M附近波動.其中示蹤粒子沿流向運動平均位移出現(xiàn)波動且波動隨時間增大的原因是：流場隨著時間的演化，本應(yīng)該聚集在中心區(qū)域和對角區(qū)域的流體微元將逐漸隨機(jī)逃逸出指定區(qū)域(具體效果可見圖16)，這些逃逸出去的粒子則反映了流場的局部不均勻性.另外結(jié)合圖3和圖6的分析結(jié)果可知，分形流場在0≤Ttrack≤100Tr內(nèi)，大尺度結(jié)構(gòu)逐漸移動到流場的邊界區(qū)域，這就導(dǎo)致了該區(qū)域流體粒子沿流向運動的擾動強(qiáng)度比Run1的中心區(qū)域更加劇烈.說明在該統(tǒng)計時間段內(nèi)，分形流場依然受到局部大尺度運動的影響，所以圖15(a)中FGT的結(jié)果曲線明顯比RGT波動的幅度更大.不過可以確定的是，當(dāng)FGT和RGT演化足夠長的時間后，兩類流場都將逐漸變成均勻各向同性湍流HIT(Homogeneous Isotropic Turbulence)狀態(tài).而在這個過程中，流場逐漸變得高度均勻，示蹤粒子沿流向運動的平均位移波動也會逐漸減小，最終等于圖11所示的理論值〈X+〉10 000/M=4.

圖15(b)與圖13(b)不同的是，無論是分形流場還是規(guī)則流場中的示蹤粒子，其沿法/展向運動的平均位移從初始的6M位置處逐漸減小，意味著Run2對角區(qū)域的流體粒子在朝著流場的中心區(qū)域運動.但很明顯的是，分形流場示蹤粒子沿法/展向運動的平均位移幅度均大于規(guī)則流場.規(guī)則流場可能是因為中心區(qū)域特殊的湍流特性，導(dǎo)致了對角區(qū)域的流體粒子朝著中心區(qū)域運動.而在分形流場中，除了中心區(qū)域特殊的湍流特性外，可能還因為邊界處的較小尺度尾流在運動過程中容易被卷入大尺度尾流中，受到大尺度運動的影響，導(dǎo)致了示蹤粒子沿法/展向運動能力要比規(guī)則流場更強(qiáng).圖15(c)的粒子沿流向運動路程與圖13(c)一致，不再贅述.圖15(d)中，分形流場和規(guī)則流場的示蹤粒子沿法/展向運動路程變化趨勢與圖13(d)幾乎一樣，分形流場的示蹤粒子沿展向運動路程大約為法向的1.1倍，而規(guī)則流場的該值約為1.2倍.分形流場中示蹤粒子沿法/展向運動路程同樣約為規(guī)則流場的3倍.

圖16 不同統(tǒng)計時刻下，Run2的分形流場和規(guī)則流場中某Y-Z平面上示蹤粒子矢量分布情況

圖16表示不同統(tǒng)計時刻下，分形流場和規(guī)則流場中某Y-Z平面上流體微元運動的二維矢量圖.其中Y-Z平面的定義為：當(dāng)示蹤粒子沿流向運動位移值在范圍3.9M≤X+≤4.1M內(nèi)，認(rèn)為示蹤粒子出現(xiàn)在該Y-Z平面上.統(tǒng)計時刻Ttrack=0時，對比圖16(a)和圖16(c)，可以看出示蹤粒子集中分布在對角區(qū)域，但分形流場中流體微元分布更為混亂.然而當(dāng)統(tǒng)計時間為Ttrack=100Tr時，盡管圖16(d)的規(guī)則流場中示蹤粒子分布展現(xiàn)出一定程度上的混亂，但沿法/展向運動還是主要集中在范圍4M≤Y+，Z+≤8M內(nèi).而此時圖16(b)的分形流場中，示蹤粒子沿法/展向運動明顯更加混亂，其中在范圍4M≤Y+，Z+≤8M內(nèi)示蹤粒子運動出現(xiàn)了明顯的漩渦狀，該位置正好對應(yīng)中心區(qū)域的第一階大尺度分形結(jié)構(gòu).同時可以看出，分形流場有大量示蹤粒子逃逸出對角區(qū)域朝著中心區(qū)域運動，結(jié)合圖15(b)的分析可知，分形流場中心區(qū)域特殊的湍流特性會對對角區(qū)域的流體微元運動產(chǎn)生影響，隨著時間演化，邊界上的小尺度尾流將逐漸被卷入大尺度尾流中，從而增大流場湍流混合度.

通過標(biāo)記流場中足夠數(shù)量的流體粒子后，對示蹤粒子的運動軌跡進(jìn)行統(tǒng)計發(fā)現(xiàn)，分形流場的對角區(qū)域主要集中為較小尺度運動，流體運動容易受到中心區(qū)域大尺度的影響，導(dǎo)致流體微元朝著中心區(qū)域運動的幅度更大.同時由于受到局部大尺度運動影響，分形流場湍流混合能力更強(qiáng)，示蹤粒子沿法/展向擴(kuò)散的程度及范圍都比規(guī)則流場更大.此外，也說明分形結(jié)構(gòu)能增強(qiáng)流場的橫向湍流混合擴(kuò)散強(qiáng)度，可以基于此原理，利用分形結(jié)構(gòu)設(shè)計，開發(fā)合理的工業(yè)元件，以實現(xiàn)對流動的被動控制.

4 結(jié) 論

本文針對具有相同阻塞率的FGT和RGT，當(dāng)流場演化隨時間至近似均勻狀態(tài)后，通過在兩類流場的不同區(qū)域釋放被動標(biāo)記的流體粒子，對示蹤粒子的運動規(guī)律進(jìn)行統(tǒng)計，具體分為兩種工況Run1和Run2.主要結(jié)論如下：

(1)當(dāng)示蹤粒子初始在某方向的釋放范圍為[0，8M]或[2M，6M]時，其粒子沿該方向運動的平均位移值將等于理論值〈X+〉10 000/M=4.另外，示蹤粒子沿流向運動平均路程與統(tǒng)計時間同樣滿足理論計算的線性關(guān)系〈SX〉i≈0.36ΔTtrack.

(2)工況Run1表明，在中心區(qū)域兩類流場的示蹤粒子沿法/展向運動的平均位移均在4M附近波動，但分形流場的局部不均勻性導(dǎo)致波動幅度要略大于規(guī)則流場.此外，示蹤粒子沿法/展向運動的平均路程大約是規(guī)則流場的3倍.

(3)工況Run2表明，分形流場的對角區(qū)域存在大尺度運動，導(dǎo)致示蹤粒子沿流向運動位移的波動幅度要大于規(guī)則流場.同時在兩類流場中，對角區(qū)域的示蹤粒子沿法/展向均朝著中心區(qū)域運動，但分形流場的粒子運動范圍更大.由于分形流場依然受到大尺度運動的影響，使得示蹤粒子沿法/展向運動的平均路程同樣滿足規(guī)則流場的3倍關(guān)系.