慣性顆粒和壁湍流的相互作用：湍流調(diào)制及顆粒分布

王貴全

(荷蘭特文特大學 流體物理課題組，恩斯赫德 7500AE，荷蘭)

0 引 言

慣性粒子在壁面湍流中的流動可用于研究沙粒在河流中的沉積、雨雪的沉積速度與分布、沙塵在地球以及火星近地面層中的傳輸以及空氣污染物在城市上空的分布等很多與地球物理、工程領(lǐng)域以及環(huán)境科學相關(guān)的重要問題[1-4]。實驗研究主要側(cè)重于流體/粒子的平均速度和雷諾應(yīng)力，以及粒子的分布，然而由于顆粒聚集與群聚效應(yīng)以及壁面的影響，導致實驗研究的難度增加，限制了實驗研究的發(fā)展[5]。隨著計算能力的不斷提高，直接數(shù)值模擬壁面湍流耦合慣性點粒子可以達到中等湍流雷諾數(shù)Reτ～O(103)（根據(jù)壁面黏性尺度），耦合顆粒數(shù)目可以達到O(108)以上，可以看到直接數(shù)值模擬已經(jīng)成為一個有效的技術(shù)手段被應(yīng)用于探究湍流粒子兩相流動這個非常復(fù)雜的物理問題，其不僅可以捕捉到最小的流動尺度，分析顆粒在湍流中的分布狀態(tài)，也可以用于研究粒子對湍流的調(diào)制。

對于慣性粒子，本文假設(shè)顆粒密度遠大于流體相密度，顆粒尺度遠小于湍流最小耗散尺度，在上述假設(shè)條件下，慣性粒子可被近似為點粒子，這樣可以極大地簡化顆粒相的數(shù)值模型。而隨著顆粒相濃度的下降，粒子與粒子間的相互作用力以及粒子對流體間的反作用力也可以逐漸的被忽略，這樣同樣也會提高計算效率。Elghobashi[6]根據(jù)顆粒相濃度的不同，將粒子與流動的耦合分為：四向耦合（同時考慮粒子對流體的反饋力以及粒子之間的相互作用力）、雙向耦合（不考慮粒子之間的相互作用力）以及單向耦合（不考慮粒子對流體的反饋力以及粒子之間的相互作用力）。本文應(yīng)用直接數(shù)值模擬和拉格朗日點粒子跟蹤方法對慣性粒子在壁面湍流這一復(fù)雜問題進行研究，主要關(guān)注粒子對湍流的調(diào)節(jié)機制及粒子分布。對湍流的調(diào)制主要研究粒子慣性對湍流內(nèi)層大尺度結(jié)構(gòu)LSMs以及外層超大尺度結(jié)構(gòu)VLSMs的影響。對粒子分布的研究包含兩部分：一是粒子聚集，這里指粒子在壁面垂直方向由于湍流場導致的濃度不均勻；二是粒子群聚形態(tài)，這里指粒子在水平截面上由于湍流場導致的不均勻分布形態(tài)。本文是對作者以往開展的相關(guān)研究工作進行小結(jié)[7-13]，結(jié)構(gòu)安排如下：第1節(jié)介紹此工作采用的數(shù)值方法及驗證過程（單向、雙向及四向耦合），第2節(jié)討論粒子對湍流內(nèi)層LSMs與外層VLSMs的調(diào)制機理（雙向耦合），第3節(jié)分析粒子在LSMs與VLSMs中的分布（雙向耦合），第4節(jié)關(guān)注粒子在重力驅(qū)動下的沉積速度（單向耦合）。

1 數(shù)值模擬方法及驗證

1.1 數(shù)值方法

1.1.1 流體相

采用直接數(shù)值模擬對不可壓縮流體求解，流向（x）和展向（z）采用周期性邊界條件并應(yīng)用譜方法求解，壁面垂直方向（y）采用二階有限差分格式，時間推進采用三階Runge-Kutta格式。不可壓縮條件通過求解壓力泊松方程來滿足。求解的連續(xù)方程和動量方程分別為：

其中ui為流體速度，p為壓力，F(xiàn)i為粒子對流體相的反作用力，v是 流體動力學黏度， ρf是流體的密度。

1.1.2 顆粒相

當粒子的密度遠大于流體相密度、粒子的尺度遠小于最小的湍流耗散尺度，粒子被近似為點粒子。此時只考慮Schiller-Naumann曳力模型，粒子的速度和軌跡可以通過積分如下方程：

在研究壁面湍流耦合慣性粒子這一問題時，最重要的控制參數(shù)為粒子的斯托克斯數(shù)（Stokes number,St），定義為粒子的特征時間尺度（ τp）和流場特征時間尺度（ τf）的比值。而由于壁面湍流邊界層內(nèi)部區(qū)域與外部區(qū)域存在不同的流場特征時間尺度，為了研究粒子與不同湍流尺度的耦合作用，定義了如下三個粒子斯托克斯數(shù)：

1.2 模型驗證[7]

在二維耦合模型中，點粒子方法是將點力耦合到最近的歐拉網(wǎng)格點上，忽略粒子體積對流體的擾動，除了上述傳統(tǒng)的點粒子方法，Capecelatro[14]介紹了一種體積平均方法，在計算粒子力的時候，不僅考慮曳力模型，同時考慮粒子大小對本地壓力和黏性應(yīng)力的影響，通過高斯核函數(shù)耦合到歐拉網(wǎng)格點上。Gualtieri[15]介紹了一種正則化點粒子方法，在計算中通過多極展開獲得球的斯托克斯流動的一階偶極子來近似點粒子對遠場的擾動。上述兩種或其他的改進點粒子方法一般都是嘗試通過引入模型來考慮粒子的體積效應(yīng)，但弊端是大大增加了計算量。為了驗證傳統(tǒng)的點粒子方法在模擬壁面湍流耦合粒子問題時的準確性，我們對傳統(tǒng)點粒子方法與體積平均方法[14]以及實驗結(jié)果[16]進行了比較，驗證過程及參數(shù)如下：

（1） 首先與已發(fā)表的采用傳統(tǒng)點粒子方法的結(jié)果[17]進行比較，目的是驗證代碼編寫的正確性。參數(shù)為Reτ= 180、St+= 30、雙向耦合。與文獻[17]比較了流體相與粒子相的平均速度及脈動速度，得到了幾乎一致的結(jié)果；

（2） 進一步比較傳統(tǒng)點粒子方法與體積平均方法[14]，來了解忽略顆粒體積效應(yīng)對模擬結(jié)果的影響。參數(shù)為Reτ= 630、St+= 2030以及Reτ= 227、St+= 58.6，四向耦合。與文獻[14]比較了流體相與粒子相的平均速度及脈動速度，得到了幾乎一致的結(jié)果；

（3） 將傳統(tǒng)點粒子方法、體積平均方法[14]與實驗結(jié)果[16]進行一對一的比較。參數(shù)為Reτ= 227、St+=58.6，模擬中采用四向耦合。比較了粒子濃度曲線，粒子聚集與群聚等穩(wěn)態(tài)值（箱形計數(shù)法，馮洛諾伊圖，概率密度函數(shù)，徑向分布函數(shù)），發(fā)現(xiàn)兩種模擬得到的結(jié)果近乎一致。但與實驗結(jié)果相比，尤其是在近壁面處，粒子的壁面垂直脈動速度以及粒子的濃度有差別，可能的原因有—數(shù)值模擬中沒有正確考慮湍流湍泳以及近壁面升力的影響[18]，或者實驗中靜電效應(yīng)以及近壁面由于粒子聚集導致測量不夠精確。對點粒子模型在近壁面處的修正是個一直沒有解決的問題，一種可能適用的方法是將小粒子看作有限大小粒子進行捕捉界面研究，但是小粒子的尺度遠小于湍流最小耗散尺度，這樣做會導致難以承受的計算網(wǎng)格數(shù)量，據(jù)作者所知，目前還缺乏相關(guān)的研究工作。

綜上，在粒子壁面湍流中，在St+= O(10～103)區(qū)間內(nèi)，傳統(tǒng)的點粒子方法和改進的點粒子方法獲得的粒子的速度與濃度一階、二階穩(wěn)態(tài)量以及流體相的速度一階、二階和三階穩(wěn)態(tài)量幾乎一致。由于傳統(tǒng)點粒子方法可以節(jié)省大量的計算時間，使得在直接數(shù)值模擬中可以引入更多數(shù)量的粒子，因此在本文的研究中，將采用這種數(shù)值方法。

1.3 計算域尺度選擇[8]

湍流內(nèi)層的大尺度結(jié)構(gòu)（LSMs）和湍流外層的超大尺度結(jié)構(gòu)（VLSMs）的存在導致流體的一階與二階穩(wěn)態(tài)值依賴于計算域的尺寸[19]。圖1展示了瞬時流向脈動速度的云圖在湍流內(nèi)層與外層、四種計算域尺寸下的模擬結(jié)果：超大尺寸Lx= 12πh，Lz= 3πh；大尺寸Lx= 6πh，Lz= 3πh；中尺寸Lx= 6h，Lz= 3h；小尺寸Lx= 2.5h，Lz= 1.5h。在湍流內(nèi)層（y+= 55處水平截面），可以看到，小計算域圖1（d）只能捕捉一對高低速條帶（高低速條帶是由于流向渦結(jié)構(gòu)將近壁面低速流體“上揚”到遠壁面的位置，或?qū)⑦h壁面高速流體“下掃”到近壁面的位置），中等以上的計算域圖1（a-c）可以得到大于三對以上的高低速條帶。在湍流外層（y+= 275處水平截面）：小計算域得不到正確的高低速條帶；中等計算域在展向方向可以捕捉到一對高低速條帶，但是流向方向有截斷；大尺度計算域的高低速條帶幾乎與超大計算域的結(jié)果相同，在流向方向幾乎沒有截斷。

在加入粒子后，粒子一階（平均速度/平均濃度）、二階穩(wěn)態(tài)（脈動速度）也依賴于計算域的尺寸[20]，在已發(fā)表的研究成果中，所采用的計算域尺寸也不盡相同（文獻[8]的表1列出了文獻中常用的計算域尺寸），可以看到，研究計算域尺寸對粒子一階與二階穩(wěn)態(tài)值的影響是非常重要的。我們對Reτ= 550、Reτ= 950以及St+= 2.41～908進行了模擬比較[8]。結(jié)果顯示在單相流中，中等尺寸計算域、大尺寸計算域和超大尺寸計算域的流體的一階、二階穩(wěn)態(tài)值幾乎相同，而小尺寸計算域會低估流體的一階、二階穩(wěn)態(tài)值。當加入粒子以后，流體相二階穩(wěn)態(tài)量在小尺寸計算域與大尺寸計算域相差較大，顆粒相的近壁面濃度在中尺寸計算域會被低估（10%～20%），尤其是對中高慣性的粒子（St+= 24.2～182），但顆粒相的二階穩(wěn)態(tài)在中尺寸計算域和大尺寸計算域得到的結(jié)果幾乎一致。因此，想要精確得到粒子在湍流內(nèi)層的濃度，需要采用大尺寸計算域。而對于模擬粒子在湍流外層的濃度，中等尺度計算域能得到和大尺度計算域相似的結(jié)果。

圖1 兩個壁面垂直位置的橫截面的瞬時流向脈動速度云圖，y+ = 55和y+ = 275分別代表湍流內(nèi)層與外層，（a-d）對應(yīng)于超大尺寸，大尺寸，中尺寸以及小尺寸[8]Fig. 1 Contours of instantaneous streamwise velocity fluctuation at two xz planes (y+ = 55 and y+ = 275) representing the inner and outer layers. Panels (a-d) refer to single-phase flow with four decreasing domain sizes, respectively [8]

2 慣性粒子對壁面湍流的調(diào)制機理

當慣性粒子的濃度達到可以對流體相產(chǎn)生影響的時候，被調(diào)制的流體相將會反之影響粒子的分布與聚集[5]。在壁面湍流中，LSMs和VLSMs分別控制湍流內(nèi)層和外層的流體動力學過程，因此我們將關(guān)注粒子對LSMs和VLSMs的調(diào)制機理[9-10]。為了突出粒子慣性的影響，在本節(jié)的計算中沒有考慮重力的影響。

2.1 對LSMs的調(diào)制機理[9]

一般地，湍流調(diào)制依賴于粒子大小和湍流尺度的比值，同時依賴于粒子的響應(yīng)時間尺度和流動特征時間尺度的比值。Gore[21]提出了用粒子與湍流積分尺度比來判斷湍流抑制還是增強。Tanaka[22]提出了用流動“含能尺度”和粒子斯托克斯數(shù)構(gòu)造的粒子動量數(shù)來預(yù)測湍流調(diào)制。盡管沒有統(tǒng)一的機理來解釋粒子的湍流調(diào)制過程，但在湍流內(nèi)層中，小慣性粒子傾向于促進湍流、大慣性粒子傾向于抑制湍流，這個推論被廣泛認可。我們通過研究粒子對近壁面擬續(xù)結(jié)構(gòu)的調(diào)制作用，對其動力學演化過程進行細致分析，目標是嘗試解釋粒子對湍流的調(diào)制機理[9]。

在湍流內(nèi)層，近壁面相干結(jié)構(gòu)的演化是一個復(fù)雜的非線性過程，對其機理的研究一直持續(xù)到現(xiàn)在。其中Hamilton[23]提出的湍流自維持過程被廣泛的用于研究LSMs的演化過程。Hamilton[23]應(yīng)用了Jiménez[24]提出的最小湍流單元概念到平板Couette流動中，這樣做的好處是可以控制展向計算域來孤立出一組大尺度條紋結(jié)構(gòu)（large-scale streaks，LSS）和大尺度渦結(jié)構(gòu)（large-scale vortex，LSV），這樣就可以把注意力放在處于同一組自維持過程中的LSS和LSV的相互轉(zhuǎn)換過程上，而不受處于其他組自維持過程中的LSS和LSV的影響，同時方便做能譜分析。而采用平板Couette流動來研究湍流自維持過程卻不采用Poiseuille流動的好處是，在上下壁面之間只存在一組湍流自維持過程，即只存在一組LSS和LSV[28]。

圖2展示了湍流自維持過程發(fā)生的位置及三種結(jié)構(gòu)的演化過程。自維持過程發(fā)生于湍流內(nèi)層中，包括三種結(jié)構(gòu)：大尺度渦結(jié)構(gòu)（LSV），大尺度條帶結(jié)構(gòu)（LSS），以及流向相關(guān)條帶結(jié)構(gòu)（x-dependent streaks）。三種結(jié)構(gòu)之間的三個轉(zhuǎn)換過程分別為：LSS形成，LSS分解，LSV再生。假設(shè)從LSV（流向不相關(guān)）開始，流體由于LSV的作用將流體“上揚”或“下掃”，形成了本地高速/低速的條紋結(jié)構(gòu)LSS（流向不相關(guān)），而LSS是不穩(wěn)定的，其失穩(wěn)產(chǎn)生流向相關(guān)的條帶，流向相關(guān)條帶引入流向速度梯度（?ux/?x），其與流向渦量（ωx）相互作用，對應(yīng)于渦量方程中的渦旋拉伸項，渦旋拉伸項的增加會促進LSV的生成，進而形成了一個完整的湍流自維持過程。應(yīng)該說明對湍流自維持過程的研究是在譜空間中進行的，上述簡單的物理描述相對應(yīng)的數(shù)學表達式請參考文獻[9,23]。

圖2 壁面湍流的內(nèi)外分區(qū)[25-26]以及自維持過程[23,27-28]. 在湍流Couette流動的最小單元中，右上的云圖顯示了“上揚”和“下掃”的LSS，向量圖顯示了一對旋轉(zhuǎn)方向相反的LSV[9]Fig. 2 Sketch of various wall regions and boundary layers in wall units[25-26] and the regeneration cycle sub-steps[23,27-28]. In the top-right is the flow field in turbulent pCf(plane Couette flow) in miniunit. Large-scale vortex (LSV) is shown by vector fields and colored by scalar value of the streamwise vorticity. Large-scale streak (LSS) is shown by the contour of the streamwise velocity magnitude[9]

近三十年的研究證明湍流自維持過程可以被用來闡述近壁面LSMs的演化機理，因此我們[9]對粒子慣性及濃度對湍流自維持過程的調(diào)制進行了細致的研究，發(fā)現(xiàn)低慣性粒子（Stturb= 0.056，用LSV的周轉(zhuǎn)時間尺度對粒子響應(yīng)時間無量綱化）在LSV中的滯留時間更長，高慣性粒子（Stturb= 0.625）會比較均勻的分布在湍流場中。相對應(yīng)的，低慣性粒子會促進湍流并縮短湍流自維持過程的周期，高慣性粒子會抑制湍流并增長湍流自維持過程的周期。由于低慣性（/高慣性）粒子的存在會增加（/減少）LSV的強度，導致近壁面湍流的增強(/減弱)。

以上，我們通過分析慣性粒子對湍流自維持過程的調(diào)制機理來研究粒子對近壁面大尺度結(jié)構(gòu)的影響。然而在實際應(yīng)用中，例如風沙流流動的雷諾數(shù)可以高達Reτ= O（106）[29]，導致湍流的內(nèi)層區(qū)間（y+＜100）只占據(jù)了整個湍流空間的很小一部分（～100/Reτ），因此除了研究粒子對近壁面LSMs的調(diào)制機理，有必要研究更高雷諾數(shù)下慣性粒子對壁面湍流外層的調(diào)制過程。

2.2 對VLSMs的調(diào)制機理[10]

隨著流動雷諾數(shù)的增加，實驗中發(fā)現(xiàn)了湍流中的另外一種重要的特征結(jié)構(gòu)，即超大尺寸結(jié)構(gòu)(VLSMs)[30]。近年來由于直接數(shù)值模擬可以達到更高的雷諾數(shù)，VLSMs也成為了數(shù)值研究湍流方向的熱點[31-33]。VLSMs的流向長度可以達到20倍壁面高度而且含有很高的能量，可以攜帶40%～65%的湍動能和30%～50%的雷諾應(yīng)力[34-35]。由于高雷諾數(shù)湍流中外層占據(jù)大部分的空間，因此研究慣性粒子對VLSMs的影響有重要的應(yīng)用價值，比如沙塵對大氣近壁面層湍流的影響[2]。然而到目前為止，慣性粒子的湍流調(diào)制數(shù)值模擬研究還局限在相對小的雷諾數(shù)(Reτ＜ 300)和對湍流內(nèi)層結(jié)構(gòu)的研究，因此我們在文獻[10]中研究了中等雷諾數(shù)(Reτ= 550，Reτ= 950)湍流中粒子慣性對VLSMs的影響。

圖3展示了慣性粒子對VLSMs調(diào)制可能存在兩個路徑：1）直接路徑，慣性粒子通過影響湍動能傳輸方程直接調(diào)制VLSMs；2）間接路徑，慣性粒子通過調(diào)制湍流內(nèi)層LSMs，進而調(diào)制湍流外層VLSMs。在譜空間里，直接路徑代表慣性粒子直接調(diào)制與VLSMs同尺度的湍流結(jié)構(gòu)，可以通過譜方法來分析湍動能傳輸方程進行研究；間接路徑代表慣性粒子通過調(diào)制內(nèi)層LSMs來達到調(diào)制外層VLSMs的目的。然而對于間接路徑，目前LSMs和VLSMs的相互作用還尚有爭論。一個觀點是LSMs和VLSMs的形成是由不同的機理導致的，理由是Jiménez[25]發(fā)現(xiàn)了LSMs的自維持機理不需要VLSMs的參與，Hwang[31]通過大渦模擬的概念進行濾波測試，發(fā)現(xiàn)VLSMs可以在沒有LSMs的存在時產(chǎn)生，Rawat[32]證明VLSMs也是一種自維持結(jié)構(gòu)并且不從LSMs吸收能量；另一種觀點是VLSMs和LSMs相互耦合存在，理由是文獻[30,36]認為VLSMs并不是一種新的結(jié)構(gòu)，而是由LSMs發(fā)展到湍流外層疊和而成，Toh[37]發(fā)現(xiàn)LSMs和VLSMs相互作用并彼此促進的周期過程，Lee[33]在湍流內(nèi)層發(fā)現(xiàn)了能量反串的過程，即LSMs向VLSMs傳遞能量。

圖3 慣性粒子調(diào)制VLSMs的兩個路徑，云圖顯示流向截面的瞬時速度，對應(yīng)Reτ = 550的明渠流動[10]Fig. 3 Schematic of two routes of inertial particle effects on VLSMs through direct impact on turbulent kinetic energy transportation or via indirect upscale energy transfer from LSMs. On the right are streamwise velocity contours in the cross-stream plane of Reτ = 550 open channel flow and the flow region[10]

我們研究[10]發(fā)現(xiàn)VLSMs增加與粒子慣性之間的非單調(diào)變化關(guān)系，即小慣性粒子（St+= 2.42根據(jù)湍流內(nèi)層時間尺度）和大慣性粒子（Stout= 6.0，8.2，根據(jù)湍流外層時間尺度）會增強VLSMs的強度（VLSMs強度通過計算預(yù)混譜中的雙模態(tài)得到），然而中等慣性粒子對VLSMs的能譜幾乎沒有影響。為了驗證上述提出的雙路徑調(diào)制VLSMs的假設(shè)，首先對小慣性粒子與大慣性粒子與流動的耦合進行人為控制，即：

（1）小慣性粒子與湍流內(nèi)層雙方向耦合，但與湍流外層單方向耦合；

（2）小慣性粒子與湍流內(nèi)層單方向耦合，但與湍流外層雙方向耦合；

（3）大慣性粒子與湍流內(nèi)層雙方向耦合，但與湍流外層單方向耦合；

（4）大慣性粒子與湍流內(nèi)層單方向耦合，但與湍流外層雙方向耦合。

測試結(jié)果為僅有（1）與（4）會增加VLSMs的強度。此測試說明，小慣性粒子促進VLSMs是通過其在湍流內(nèi)層的反饋力，而大慣性粒子促進VLSMs是通過其在湍流外層的反饋力。進一步，用譜方法分析了雷諾應(yīng)力輸運方程，發(fā)現(xiàn)高慣性粒子對流體的反饋力作用于雷諾剪切應(yīng)力輸運方程，促進與VLSMs相同尺度的雷諾剪切應(yīng)力的生成，雷諾剪切應(yīng)力作用于流向湍動能輸運方程而直接增加VLSMs的強度。物理解釋為：VLSMs與高慣性粒子相互作用，粒子的反饋力促進與VLSMs尺度相同的“上揚”和“下掃”流動，隨后這些超大尺度的“上揚”和“下掃”結(jié)構(gòu)從平均流動中吸收能量，最終直接促進同等尺度VLSMs的生成；低慣性粒子增強湍流自維持過程[9]并增加LSMs的強度，通過能量反向級串，最終促進VLSMs的生成。這種間接調(diào)制機理，即小粒子通過調(diào)制LSMs并最終調(diào)制VLSMs，支持了上述的第二種觀點[30,33,36-37]。

以上討論了慣性粒子的存在對LSMs和VLSMs的調(diào)制機理，除了粒子對流場的影響以外，研究粒子在LSMs和VLSMs的聚集狀態(tài)也對大渦模擬粒子流的準確性及局限性提供借鑒意義。由于計算能力的限制，目前直接數(shù)值模擬耦合大量的慣性粒子（O（108））只能達到中等的雷諾數(shù)（Reτ＜1 000），遠遠不能達到實際情況中的超高雷諾數(shù)Reτ～O（106）。在高雷諾數(shù)的情況下，由于大渦模擬過濾掉小于計算網(wǎng)格的流動結(jié)構(gòu)，導致計算中得到的粒子本地處的流體速度并不精確（即滑移速度不準確），Marchioli[38]發(fā)現(xiàn)即使是在低雷諾數(shù)Reτ= 180的情況下，大渦模擬也會低估粒子近壁面的聚集，同時粒子的群聚形態(tài)也不同于直接數(shù)值模擬。而在高雷諾數(shù)的情況下，采用大渦模擬往往只能捕捉到VLSMs而過濾掉LSMs，想要準確地應(yīng)用大渦模擬到高雷諾數(shù)湍流粒子流中[39-40]，有必要衡量這兩種重要流動特征結(jié)構(gòu)對不同慣性的粒子分布的影響。

3 粒子在LSMs及VLSMs中的聚集與

采用了如下的策略來應(yīng)用直接數(shù)值模擬分離LSMs與VLSMs對粒子的影響：

（1）應(yīng)用截斷計算空間的方法分離出LSMs，比較粒子分布在只有LSMs及在同時有LSMs和VLSMs時的不同，見圖1；

（2）在計算過程中過程中，在譜空間分離出LSMs和VLSMs，人為的控制粒子與這兩種結(jié)構(gòu)流場的耦合，比較粒子分布在只有VLSMs及同時有LSMs和VLSMs時的不同，見圖4。

圖4 Reτ = 550單相流，y+ = 100水平截面以及四個邊界截面的瞬時流向脈動速度云圖，“+”代表依據(jù)黏性尺度進行無量綱化[11]Fig. 4 Instantaneous contours of streamwise velocity fluctuation on a wall-parallel plane at y+ = 100(and domain boundary walls) in single-phase flow，normalized by uτ[11]

對于測試（1），在第1.3節(jié)中已經(jīng)介紹了四種計算域尺寸對粒子的影響，重點關(guān)注粒子的一階與二階穩(wěn)態(tài)[8]。本節(jié)用截斷計算域作為一種分離LSMs的方法，關(guān)注VLSMs的存在對粒子分布的影響。首先根據(jù)文獻[41]，選擇計算域流向尺寸足夠大到包含一個LSMs并足夠小到排除外層VLSMs對內(nèi)層LSMs的影響，在Reτ= 550時，我們選擇截斷計算域尺寸為Lx+= 1 375、Lz+= 825，而大尺寸計算域尺寸為Lx+=10 367、Lz+= 3 456，通過比較二維空間的預(yù)混譜，截斷計算域得到的湍流內(nèi)層能量譜與大尺寸計算域以及文獻[42]的結(jié)果幾乎吻合，但湍流外層的能量譜由于截斷域的影響，小尺寸計算域由于尺寸太小顯然捕捉不到能量譜的雙模態(tài)形式。在加入慣性粒子（低慣性St+= 24.2和高慣性St+= 182）以后：一方面，在湍流內(nèi)層尤其是接近壁面處，粒子的濃度小于大尺寸計算域接近20%，這歸因于外層VLSMs的存在增強了湍流湍泳（由于脈動速度在壁面垂直方向梯度誘導的粒子運動，粒子趨向于向湍流度弱的方向移動），由于粒子的平均濃度在兩種計算域中保持相同，因此截斷域的湍流外層粒子濃度則小于大計算域中的粒子濃度；另一方面，我們通過馮洛諾伊圖來計算粒子在壁面垂直方向上的群聚形態(tài)[43]，發(fā)現(xiàn)無論對低或高慣性粒子，截斷域總是趨向于低估粒子的群聚效應(yīng)，尤其是在湍流外層區(qū)域（y+＜ 100），這同樣歸因于外層VLSMs的影響。

對于測試（2），在計算過程中我們?nèi)藶榈目刂屏Ｗ优c流場特定結(jié)構(gòu)的耦合，具體的步驟為：將速度場變換到譜空間，將流場分為LSMs和VLSMs兩部分，再從譜空間變換到速度場（見圖4），粒子只與LSMs或VLSMs對應(yīng)的速度場耦合，這種人為控制的粒子-流場耦合在計算的每一步都進行一次，因此計算時間較長。通過以上方式可以得到三種粒子分布，即全流場下、LSMs流場下、VLSMs流場下的粒子分布（見圖5），可以看到對于小慣性粒子，其粒子群聚主要發(fā)生在湍流內(nèi)層，直觀上LSMs流場下粒子湍流內(nèi)層的群聚看起來與全模擬相似（圖5（a）與（c）比較）。而在湍流外層，由于其響應(yīng)時間尺度遠小于特征結(jié)構(gòu)VLSMs的時間尺度，因此小慣性粒子在湍流外層中趨向于無慣性粒子；對于高慣性粒子，由于其響應(yīng)時間尺度與VLSMs的特征尺度相當，因此高慣性粒子的群聚發(fā)生于湍流外層，然而即使是直觀上，VLSMs流場下的粒子在湍流外層的群聚形態(tài)也和全模擬相差很大（圖5（b）與（d）比較）。隨后對粒子聚集（粒子濃度曲線）以及粒子群聚（馮洛諾伊圖計算）進行了比較，發(fā)現(xiàn)在湍流內(nèi)層，LSMs流場會低估慣性粒子的濃度，在湍流外層，VLSMs流場可以得到和全流場相近的慣性粒子的濃度。對于粒子群聚現(xiàn)象，LSMs流場結(jié)果更接近于全流場結(jié)果，而VLSMs流場結(jié)果與全流場結(jié)果相差很大。

根據(jù)測試（1）與測試（2），發(fā)現(xiàn)LSMs影響粒子的聚集并控制粒子的群聚。在大渦模擬粒子流中，過濾掉LSMs可以估計湍流外層的粒子濃度，但會低估湍流內(nèi)層的粒子濃度。大渦模擬捕捉不到正確的粒子群聚，導致粒子與流場的雙向耦合作用到不合適的流場結(jié)構(gòu)中，因此不適宜采用大渦模擬和粒子的雙向耦合研究粒子的反饋作用對流場的調(diào)制。

圖5 粒子在三種流場中的分布[11]Fig. 5 Instantaneous snapshots of particle locations (black dots)[11]

4 影響粒子沉積的機理研究[12-13]

在一些地球物理相關(guān)測量中，能準確估算雪、雨、沙塵等的表面沉積速率對研究降水量、土地沙化、城市污染有重要的實際應(yīng)用價值[44-45]。這些粒子的沉積過程，可以看作慣性粒子在重力作用下穿越壁面湍流場到達壁面這個簡化的物理過程。當粒子的St+～0時，粒子可以近似看做具有沉降速度的被動標量，Rouse[46]推導出了用于預(yù)測粒子垂直濃度的冪定律，冪指數(shù)與慣性粒子的斯托克斯沉積速度成一定比例，被廣泛應(yīng)用于地球物理的相關(guān)研究中。然而對于慣性粒子，文獻[46]給出的冪定律則不再適用，Richter[47]嘗試加入粒子對流速度的一階展開（見文獻[48]）去修改冪定律，通過與直接數(shù)值模擬結(jié)果比較，發(fā)現(xiàn)這樣考慮粒子慣性效應(yīng)的一階修正只使用于St+＜＜1的情況，對于St+＞1的慣性粒子，目前存在的理論模型仍然很難準確預(yù)測粒子的垂直濃度。對于慣性粒子，文獻[46]結(jié)論不適用的原因是其假設(shè)粒子在湍流場中均勻分布、粒子的沉積速度等于其斯托克斯沉積速度，這樣的假設(shè)對于St+＞1的慣性粒子不再適用，因為即使在均勻一致湍流中，由于慣性粒子趨向于遠離渦的中心，粒子傾向于被渦旋轉(zhuǎn)向重力方向加速，導致慣性粒子的沉降速度大于其斯托克斯沉積速度[49-50]。而在壁面湍流中，由于多種湍流尺度的存在，不同慣性的粒子會偏向聚集于湍流場的不同結(jié)構(gòu)中，影響粒子加速的因素變得更加復(fù)雜。

在文獻[12-13]中，我們首先通過推導相空間下的概率密度分布的主方程（phase-space master PDF equation）得到了控制粒子穩(wěn)態(tài)的輸運方程，進而推導出控制粒子沉降速度的方程（數(shù)學推導過程請參考文獻[12,51]）：

其中vp(t) 為粒子的沉降速度，y為壁面垂直方向，ρ為邊〈緣概率密度函〉 數(shù)用來描述粒子的空間分布，S=(vp(t)-〈vp(t)〉)2為粒子脈動速度的協(xié)方差，up(t)為粒子處的流體速度。通過對此控制方程的分析，發(fā)現(xiàn)粒子的沉降速度是由如下五個因素來控制，對應(yīng)方程右側(cè)由左至右依次為：粒子加速項，粒子擴散項，湍流湍泳項，粒子本地流場速度項以及重力加速項。對上述五個因素的物理描述為：粒子加速項是由粒子在壁面垂直方向的速度梯度導致的，粒子擴散項是由粒子在壁面垂直方向的濃度梯度導致的，粒子本地流場速度項是由粒子在湍流中的偏向聚集現(xiàn)象導致的，重力加速項是由外力導致的。隨后通過直接數(shù)值模擬計算了上述五個因素在零粒子通量以及有限大小粒子通量兩種情況下對粒子沉降速度的貢獻。進一步的，我們討論了擴展文獻[46]中的模型從St+～0粒子到St+＞ 1粒子的可能性，發(fā)現(xiàn)模型中的不封閉項包括由于粒子聚集導致的漂移系數(shù)以及粒子壁面垂直方向的湍動能。我們目前仍然尚無法提出滿意的模型用來預(yù)測St+＞ 1慣性粒子的壁面垂直方向的濃度曲線，這將是未來要開展的工作之一。

5 結(jié) 論

本文應(yīng)用直接數(shù)值模擬耦合拉格朗日點粒子這種相對簡單的數(shù)值模型，對慣性粒子和壁面湍流這個復(fù)雜且具有重要實際應(yīng)用的問題開展了一些相關(guān)研究。重點關(guān)注的物理問題是粒子對壁面湍流的調(diào)制，以及粒子在兩種特征湍流結(jié)構(gòu)中的分布以及粒子的沉積速度。在未來的研究中，將繼續(xù)嘗試建模預(yù)測粒子垂向濃度曲線，熱對流效應(yīng)對壁面湍流及粒子分布的影響，粒子在復(fù)雜流場結(jié)構(gòu)中的運動（例如旋轉(zhuǎn)湍流），帶電/磁粒子在壁面湍流中的運動，以及應(yīng)用大渦模擬技術(shù)研究復(fù)雜地形下粒子分布等相關(guān)研究工作。

致謝：本文涉及到的工作是作者在美國圣母大學的研究內(nèi)容，全部已公開發(fā)表到經(jīng)過同行評審的專業(yè)期刊或尚未經(jīng)過同行評審的預(yù)印版。本工作得益于美國陸軍研究實驗室以及海軍研究辦公室的項目支持，計算資源來自于高性能計算現(xiàn)代化項目（HPCMP）以及圣母大學研究計算中心。相關(guān)的研究工作得益于幾位合作者的貢獻：美國圣母大學的Hyungwon John Park博士以及David H. Richter副教授，美國杜克大學的Andrew Bragg助理教授，美國明尼蘇達大學的Kee Onn Fong博士以及Filippo Coletti副教授（現(xiàn)工作于蘇黎世聯(lián)邦理工大學），美國密歇根大學的Jesse Capecelatro助理教授。感謝幾位優(yōu)秀學者訪問本人所在的美國圣母大學課題組時所提出的寶貴建議：美國加州大學洛杉磯分校的Marcelo Chamecki教授以及美國華盛頓大學的James J.Riley教授。感謝已發(fā)表論文和本論文匿名審稿人提出的寶貴建議。最后感謝相關(guān)專家的推薦及《空氣動力學學報》編輯部的邀請。