近壁湍流中微小非球形顆粒取向行為研究綜述

崔智文，趙立豪

(清華大學(xué) 航天航空學(xué)院 工程力學(xué)系，北京 100084)

0 引 言

復(fù)雜流動中的顆粒運動在自然環(huán)境和工業(yè)生產(chǎn)等領(lǐng)域中普遍存在，比如云層冰晶以及雨雪的形成[1]、浮游生物的運動與聚集現(xiàn)象[2]、紙漿纖維與造紙過程[3-4]、化學(xué)制藥過程[5]等。同時，顆粒的形狀通常是各向異性的。但是過去大部分的研究用等效的球形點顆粒模型對顆粒運動進行描述，開展顆粒的沉積、聚集以及與湍流相互作用等問題的研究[6-10]。然而，顆粒的非球形的形狀特征在一定程度上會對顆粒的輸運產(chǎn)生影響。與此同時，在一些特定的自然與工程問題中，顆粒的形狀特征不可忽略。因此由顆粒形狀各向異性引起的取向行為就具有十分重要的研究意義，比如造紙過程中細(xì)長的紙纖維的取向影響紙張的力學(xué)性能[3-4]；纖維增強材料中纖維的排列決定了材料的微結(jié)構(gòu)，進而影響材料的力學(xué)性能；在減阻控制中，纖維的取向行為影響纖維與湍流之間的相互作用[11]；浮游生物的取向行為會影響生物的聚集與遷移過程[12]等。

Fan 和 Ahmadi[13-14]、Zhang等[15]較早采用歐拉-拉格朗日方法對桿狀顆粒在壁湍流中的運動問題進行數(shù)值模擬。隨后，Mortensen 等[16-17]、Marchioli等[6,18-19]、Marchioli & Soldati[20]和Challabotla等[21-23]對非球形顆粒在壁湍流中的輸運過程開展了大量相關(guān)研究工作。其中，歐拉-拉格朗日方法是背景流場在歐拉的觀點下進行求解，而顆粒則通過拉格朗日追蹤的方法進行求解。壁湍流是一類在壁面附近存在強剪切與豐富的擬序結(jié)構(gòu)的湍流，比如槽道、圓管以及邊界層湍流等[24]。因此，壁湍流中的非球形顆粒兩相流研究主要關(guān)注以下幾類問題：1）非球形顆粒在壁湍流中的統(tǒng)計行為（平動、轉(zhuǎn)動以及取向）；2）非球形顆粒與湍流中的相干結(jié)構(gòu)之間的相互作用；3）非球形顆粒對湍流的減阻控制的作用。

本文主要綜述了近些年關(guān)于微小非球形顆粒在壁湍流中取向行為的研究工作進展，具體結(jié)構(gòu)如下：第1節(jié)介紹相關(guān)工作使用的理論與數(shù)值模擬方法，第2節(jié)討論非球形顆粒在槽道湍流中的取向行為，第3節(jié)介紹顆粒傾向性取向行為的機理，第4節(jié)討論顆粒的取向行為會在近壁區(qū)呈現(xiàn)“形狀敏感性”的現(xiàn)象，最后對本文的內(nèi)容進行總結(jié)和展望。

1 數(shù)值方法

前期的壁湍流顆粒數(shù)值研究主要方法為歐拉-拉格朗日耦合的直接數(shù)值模擬。這里對流場及顆粒求解方法做簡單的介紹。

1.1 流體相

直接數(shù)值模擬對不可壓縮的納維-斯托克斯方程進行求解，求解的連續(xù)方程與動量方程分別為：

其中ui為 流體速度，p為壓力，v是流體的運動學(xué)黏度，ρf是流體的密度。本文主要介紹的物理模型為兩平板間的槽道湍流，其中流向（x）和展向（y）為均勻的周期方向，且在壁面處滿足速度不可滑移條件。在流向與展向均使用偽譜方法求解，而在壁面垂直方向（z）采用二階有限差分格式。在時間推進上采用顯示二階Adams-Bashforth格式。

主要算例的計算域大多為1 2h×6h×2h，其中h為槽道高度的一半。算例的計算網(wǎng)格數(shù)為1 923，在流向與展向為均勻網(wǎng)格，而壁面法向采用雙曲正切函數(shù)形式的非均勻網(wǎng)格。算例中基于壁面摩擦速度的摩擦雷諾數(shù)Reτ=180。在本文中上標(biāo)“+”代表著變量已被黏性尺度進行無量綱化處理。

1.2 顆粒相

本文主要介紹的顆粒尺寸遠小于流體的Kolmogorov尺度，且顆粒相為稀疏懸浮，因此不考慮顆粒對流體的反作用以及顆粒與顆粒之間的碰撞作用。同時，顆粒的尺寸非常?。紤]的微粒弛豫時間與流體弛豫時間小于0.1，但顆粒與流體的密度比可以達到1 000左右），顆粒慣性與流體慣性對顆粒運動學(xué)的影響可以被忽略[25]。此時，顆?？山普J(rèn)為跟隨流體的軌跡線，即

其中up,i為顆粒速度。對于跟隨流體的顆粒，顆粒速度就近似等于當(dāng)?shù)氐牧鲌鏊俣?。本文主要考慮軸對稱橢球顆粒，因此顆粒的取向（軸對稱橢球顆粒的回轉(zhuǎn)軸方向）可以通過Jeffery方程[26]得到，即

其中pi為 顆?；剞D(zhuǎn)軸的方向（單位向量）；Oij為流體的旋轉(zhuǎn)張量，即為流體的變形率張量，即而 λ為顆粒的形狀參數(shù)，即為顆粒的回轉(zhuǎn)軸與赤道軸的長度之比。有時也使用形狀因子對顆粒形狀進行描述。不同于顆粒參數(shù)λ，形狀因子 Λ的取值范圍為 [-1,1]。 當(dāng) λ＞1時 ， 0＜Λ＜1，此時顆粒為桿狀顆粒；當(dāng) 0＜λ＜1時 ， -1＜Λ＜0，那么顆粒則為碟狀顆粒；當(dāng) λ=1時，顆粒呈球形，此時Λ=0。通過顆粒的取向方程可知，顆粒的形狀參數(shù)λ（或形狀因子Λ）主要是調(diào)節(jié)流體變形率在顆粒取向方程中的占比，進而影響顆粒的取向行為。當(dāng)形狀參數(shù) λ→∞或 0時（形狀因子 |Λ|=1），此時流體變形率的作用最強。

本文中顆粒相的時間推進格式與流體相的時間推進格式保持一致，而顆粒處的流體信息通過二階三維拉格朗日插值格式從顆粒附近的流體網(wǎng)格中插值得到。同時，顆粒與壁面之間采用完全彈性碰撞模型。取向隨機的顆粒在初始時刻被均勻地以當(dāng)?shù)亓鲌鏊俣确湃氤浞职l(fā)展的槽道湍流中，待顆粒在流場中充分發(fā)展后進行統(tǒng)計與分析。

2 顆粒在壁面附近的取向行為

在槽道湍流中，因為槽道內(nèi)部湍流非均勻且各向異性，非球形顆粒的取向與轉(zhuǎn)動行為受顆粒形狀的影響較大，尤其對于極細(xì)長的桿狀顆粒和極扁平的碟狀顆粒。而且，顆粒的行為在槽道中的不同位置處會呈現(xiàn)出不同的行為特征。

Challabotla等[22]針對非球形顆粒在槽道中的取向行為展開了研究。圖1展示了瞬時細(xì)長顆粒與扁平顆粒在近壁面處的分布，可以明顯看到顆粒具有傾向性的取向，同時這種取向分布與顆粒形狀有關(guān)。圖2展示了非球形顆?；剞D(zhuǎn)軸與流向、展向與法向方向的夾角余弦絕對值的統(tǒng)計平均值。在槽道中部，因為中部的流動狀態(tài)趨近于均勻各向同性湍流[27]，所以非球形顆粒的運動行為與其在均勻各向同性湍流中的結(jié)果基本一致，即桿狀顆粒傾向性地朝著渦量方向而碟狀顆粒垂直于渦量方向。在慣性參考系中看，顆粒的回轉(zhuǎn)軸方向與慣性參考系各個軸的夾角余弦的絕對值的統(tǒng)計值均趨于0.5，這意味著非球形顆粒在槽道中部的取向在空間上是趨于一種隨機分布的狀態(tài)。然而，從槽道的黏性底層區(qū)（z+＜5）延伸到z+≈30的緩沖區(qū)，桿狀顆粒與碟狀顆粒的取向行為呈現(xiàn)出明顯的差異，即細(xì)長的桿狀顆粒在壁面附近傾向性地朝著流向方向，而碟狀顆粒的回轉(zhuǎn)軸方向傾向性地朝著壁面法向。顆?；剞D(zhuǎn)軸的方向與流向或法向方向的夾角余弦統(tǒng)計值隨顆粒形狀參數(shù)從0.01至50的變化而單調(diào)變化。該變化趨勢與非球形顆粒在各向同性湍流中顆粒與渦量夾角隨形狀參數(shù) λ變化趨勢類似。

圖1 非球形顆粒在平行壁面平面上的瞬時分布圖（z+≈10）Fig. 1 Instantaneous distribution of non-spherical particles in a wall-parallel plane at z+≈10

圖2 非球形顆粒回轉(zhuǎn)軸與流向、展向與法向方向的夾角余弦絕對值的統(tǒng)計平均值[22]Fig. 2 Average absolute cosine values between particle symmetry axis and the streamwise,the spanwise and the wall-normal direction[22]

通過Jeffery方程[26]可知，非球形顆粒的轉(zhuǎn)動行為與取向行為是一種相互影響的關(guān)系，即顆粒取向的時間導(dǎo)數(shù)其中 ω為顆粒的轉(zhuǎn)動角速度。在槽道湍流中，壁面附近的渦量場具有較強的各向異性，并且在展向存在非常強的渦量。圖3(a)展示了顆粒與流體微團的展向平均旋轉(zhuǎn)角速度。其結(jié)果表明在近壁區(qū)非球形顆粒的展向轉(zhuǎn)動受到了抑制，其抑制效果隨形狀偏離球形的程度的增大而增強。但是隨著統(tǒng)計位置逐漸遠離壁面 （遠離緩沖區(qū)），非球形顆粒的轉(zhuǎn)動行為受顆粒形狀的影響程度會逐漸減弱。非球形顆粒在壁面的轉(zhuǎn)動受抑制也進一步解釋了為什么顆粒在壁面存在較強的傾向性的取向行為。前期工作認(rèn)為非球形顆粒在近壁區(qū)的行為（尤其在黏性底層）與顆粒在線性剪切流中的現(xiàn)象類似，并認(rèn)為顆粒的轉(zhuǎn)動行為類似Jeffery軌跡[26]，即桿狀顆?；剞D(zhuǎn)軸在流向附近（或碟狀顆粒回轉(zhuǎn)軸在法向附近）時需要停留較長的時間才會突然進行翻轉(zhuǎn)。因此，在統(tǒng)計的基礎(chǔ)上可以得到非球形顆粒在壁面附近具有傾向性的統(tǒng)計學(xué)取向行為。圖3(b-d)則進一步展示了顆粒在慣性參考系中各個方向上轉(zhuǎn)動角速度的脈動值隨顆粒形狀參數(shù)的復(fù)雜變化規(guī)律。與此同時，揭育澄等人[28-29]也研究了中等雷諾數(shù)Reτ=1000中非球形顆粒的取向行為，并發(fā)現(xiàn)顆粒的取向行為的統(tǒng)計結(jié)果基本與低雷諾數(shù)Reτ=180接近[28-29]。該結(jié)果也說明了非球形顆粒在近壁處的取向行為具有一定雷諾數(shù)無關(guān)的特征[28-29]。

圖3 顆粒轉(zhuǎn)動角速度統(tǒng)計值[22]Fig. 3 Statistics of particle angular velocities[22]

3 顆粒傾向性取向行為的機理

顆粒取向行為的機理一直以來都是關(guān)注的熱點問題。在均勻各向同性湍流的研究中，目前普遍接受的觀點主要有兩種：1）非球形顆粒的取向與流體的渦量以及變形率張量的第二特征值方向相關(guān)[30-32]；2）非球形顆粒的取向與流體的拉格朗日拉伸與壓縮方向相關(guān)[33-34]。前者主要是基于歐拉的觀點去討論顆粒與流場物理量之間的關(guān)系，其中渦量隨時間演化的拉格朗日控制方程與桿狀顆粒退化的演化方程相似，唯一的不同點在于渦量方程多出了黏性擴散項[30,32]。同時，由于渦量與變形率張量之間的關(guān)聯(lián)，顆粒的取向同樣與變形率張量的第二特征值方向具有較強的相關(guān)性。然而，顆粒與渦量相關(guān)的解釋僅在各向同性湍流或遠離壁面區(qū)域的湍流有效，但是在具有較強剪切的壁面附近并不成立。因為壁面的存在，流場在壁面會存在較強的展向渦量，但是桿狀顆粒與碟狀顆粒均不會傾向性地朝著展向[22]，所以桿狀顆粒與碟狀顆粒被觀察到垂直于渦量方向。同時，在二維流場中，渦量方向永遠垂直于流場的平面。因此，第一種解釋具有一定的局限性。對于第二種觀點，極細(xì)長的桿狀顆?？梢钥闯蔁o限小的物質(zhì)線段，而該物質(zhì)線段與流體的拉格朗日拉伸方向漸近一致[30,35-36]，當(dāng)時間足夠長時，桿狀顆粒與拉格朗日拉伸方向一致。Ni 等[34]發(fā)現(xiàn)在均勻各向同性湍流中形狀參數(shù) λ＞10的桿狀顆粒與拉格朗日拉伸方向就具有較好的一致性，而且隨形狀參數(shù) λ的變化并不明顯。為了進一步驗證其在具有各向異性的湍流中是否也同樣適用，趙立豪等[37]在槽道湍流中分析了非球形顆粒與流體拉格朗日拉伸與壓縮方向之間的相關(guān)性。

3.1 拉格朗日拉伸與壓縮方向

流體的拉格朗日拉伸與壓縮方向是指初始時刻為球形的流體微團沿著拉格朗日軌跡線發(fā)生拉伸與壓縮變形的方向[34]。因此，為了得到流體的拉格朗日拉伸與壓縮方向，首先需要沿著流體跡線積分得到對應(yīng)的變形梯度張量Fij，即[33-34]

其中Aij是 流體的速度梯度張量，而Fij是流體的變形梯度張量。變形梯度張量Fij表征了流體微團沿軌跡線的變形情況。Fij不是實對稱張量，通常采用柯西格林張量來表征流體微團的變形。本文采用左柯西格林張量進行表征，即[33-34]

通常設(shè)定初始時刻的Fij=δij（ δij為Kronecker符號），即初始時刻為單位球形張量，隨時間積分一段時間后，對M進行特征分解。其中M的最大特征值對應(yīng)的特征方向作為拉格朗日拉伸方向，記為eL1。同理，最小特征值對應(yīng)的特征方向為拉格朗日壓縮方向，記為eL3。第二特征值對應(yīng)的特征方向記為eL2。而在壁面附近，由于流動近似為線性剪切流，所以在統(tǒng)計結(jié)果中拉格朗日的拉伸方向指向流向，而壓縮方向指向壁面法向[29]。

3.2 非球形顆粒與拉格朗日拉伸與壓縮方向的相關(guān)性

圖4展示了顆粒回轉(zhuǎn)軸方向p與左柯西格林張量的主軸方向eLi的 相關(guān)隨積分時間的變化，即〈 (eLi·p)2〉。若相關(guān)值為1/3，那么顆粒相對于柯西格林張量的主軸方向是隨機分布的；而相關(guān)值為1意味著顆粒與對應(yīng)的主軸方向完全一致。圖4橫軸為相對初始時刻的積分時長[37]。圖4的結(jié)果表明了桿狀顆粒逐漸趨近于拉格朗日拉伸方向，而碟狀顆粒逐漸趨近于拉格朗日的壓縮方向。圖5為顆?；剞D(zhuǎn)軸方向p與左柯西格林張量的主軸方向eLi的 相關(guān)沿空間分布，即〈 (eLi·p)2〉。圖5選取積分時間t+=+136至t+=+144進行統(tǒng)計，圖片來源于文獻[37]。通過顆?；剞D(zhuǎn)軸的方向與柯西格林張量的三個主軸方向的夾角的空間分布圖，發(fā)現(xiàn)桿狀顆粒與拉格朗日拉伸方向以及碟狀顆粒與拉格朗日壓縮方向在近壁區(qū)域存在非常強的相關(guān)性。但是，細(xì)長的桿狀顆粒和扁平的碟狀顆粒的取向與拉格朗日的拉伸和壓縮方向還是存在比較明顯的差異，至少與Ni等[34]在均勻各向同性湍流中提出來λ＞10的形狀條件并不相符。而這種差異產(chǎn)生的原因有待進一步研究。

圖4 顆粒回轉(zhuǎn)軸的方向與左柯西格林張量的三個主軸方向的夾角隨時間演化的關(guān)系[37]Fig. 4 Time evolution of the alignment of the orientation vector p of spheroidal particles with aspect ratio relative to the three eigenvectors of left Cauchy-Green tensor[37]

圖5 顆?；剞D(zhuǎn)軸的方向與柯西格林張量的三個主軸方向的夾角的空間分布圖[37]Fig. 5 Variation of the alignment of the orientation vector p of spheroidal particles with aspect ratio relative to the three eigenvectors of left Cauchy-Green tensor with different aspect ratio[37]

4 顆粒取向行為的形狀敏感性問題

在第3節(jié)中，本文介紹到桿狀顆粒與拉格朗日拉伸方向以及碟狀顆粒與拉格朗日的壓縮方向具有較強的相關(guān)性。但是，顆粒與拉格朗日拉伸與壓縮方向的差異隨顆粒形狀參數(shù) λ變化依然明顯，這與Ni 等[34]觀察到 λ＞10的桿狀顆粒在各向同性湍流中的取向就基本與拉格朗日拉伸方向一致的結(jié)論存在差別。本章主要討論在近壁面附近這種差異是如何體現(xiàn)，及其原因。

4.1 同一軌跡線上顆粒與拉格朗日拉伸方向的差異

首先，在崔智文等[38]的工作中，他們選取一條軌跡線，并截取顆粒已經(jīng)充分發(fā)展且在壁面附近有較長時間停留但與壁面無碰撞過程的軌跡段，如圖6所示。圖6展示了同一條軌跡線上形狀參數(shù) λ=23.3的桿狀顆粒的取向的時間演化圖與拉格朗日拉伸方向的差異。圖中 φ是相對于慣性參考系的方位角，θ是俯仰角。槽道的上壁面在z+=360 ，當(dāng)2 700＜t+＜3200時，顆粒處在上壁面的粘性底層區(qū)z+＞355，圖片來自文獻[38]。顆粒的形狀參數(shù) λ=23.3， 對應(yīng)于Λ=0.9963（圖中藍色虛線），此時桿狀顆粒的 Λ參數(shù)已經(jīng)非常接近于1，其與1的差別僅為 δΛ=1-Λ=0.0037。而顆粒的形狀因子 Λ=1（圖中紅色實線）意味著顆粒的長細(xì)比無限大，并且在充分發(fā)展的顆粒場中無限長細(xì)比的桿狀顆粒的取向與流體的拉格朗日拉伸方向一致。而在差別僅為 δΛ=0.0037的情況下，人們從直觀上認(rèn)為桿狀顆粒的取向應(yīng)該與拉格朗日的拉伸方向基本一致。然而，圖6的結(jié)果表明，除了遠離壁面的區(qū)域，桿狀顆粒的取向與拉格朗日拉伸方向在近壁區(qū)卻呈現(xiàn)出迥然不同的運動行為。在非常靠近壁面的位置，拉格朗日的拉伸方向基本朝著流向方向，但桿狀顆粒卻持續(xù)地翻轉(zhuǎn)。圖6中的方位角 φ是桿狀顆粒回轉(zhuǎn)軸在x-z平 面的投影與x軸 的夾角，而俯仰角 θ是桿狀顆?；剞D(zhuǎn)軸與x-z平面的夾角。

圖6 λ=23.3的桿狀顆粒（Λ= 0.996 3，對應(yīng)藍色虛線）與拉格朗日拉伸方向（Λ = 1，對應(yīng)紅色實線）沿同一軌跡線的取向行為差異[38]Fig. 6 Angular dynamics of a slender rod (λ = 23.3, Λ = 0.996 3,blue dashed lines) and the Lagrangian stretching direction eL1(Λ = 1, red solid lines) along a trajectory[38]

4.2 桿狀顆粒在拉格朗日坐標(biāo)系中的相對分布

為進一步分析，崔智文等[38]對細(xì)長桿狀顆粒在拉格朗日坐標(biāo)系（由柯西格林張量三個主軸方向組成，即eL1、eL2和eL3)的相對取向的概率密度分布進行了研究。圖7與圖8分別展示了該分布在槽道中部（z+=180） 和槽 道 黏性底層（z+=4）中 的 情 況。其中方位角α是桿狀顆?；剞D(zhuǎn)軸在eL1-eL3平面的投影與eL1的 夾 角，而 俯 仰 角 β是 桿 狀 顆 粒 回 轉(zhuǎn) 軸 與eL1-eL3平面的夾角。結(jié)果表明不管是在槽道中部還是壁面的黏性底層，其歐拉角α與 β的概率密度分布函數(shù)都存在平臺區(qū)與冪律區(qū)這兩個明顯的區(qū)域。其中平臺區(qū)的寬度隨形狀差別 δΛ=1-Λ（即無限長細(xì)比與有限長細(xì)比顆粒之間的形狀因子的差別）的增加而變寬，而冪律區(qū)中的分布函數(shù)接近P(α)～α-2。崔智文等[38]認(rèn)為平臺區(qū)的寬度反映了桿狀顆粒取向與拉格朗日拉伸方向的差異。平臺越窄意味著桿狀顆粒在絕大部分時候基本與拉格朗日拉伸方向一致。但隨著平臺的變寬，桿狀顆粒取向與拉格朗日拉伸方向的差異逐漸增強。通過圖7與圖8的對比，不難發(fā)現(xiàn)，不管在槽道中部還是在黏性底層，歐拉角α與β的分布是類似的，但是平臺的寬度在壁面附近要明顯大于槽道中部的情況。如果將槽道中各層中不同形狀顆粒的α分布的平臺的寬度αc進行提?。▓D7與圖8中的虛線），并將其與顆粒形狀 δΛ=1-Λ 進行比較（如圖9所示），其結(jié)果表明在槽道的各個層的平臺寬度αc與 δΛ呈較強的線性相關(guān)。圖9的結(jié)果說明了圖7與圖8中發(fā)現(xiàn)的規(guī)律不僅存在于壁面附近存在強剪切的區(qū)域，而且也適用于槽道的各個區(qū)域。

圖7 槽道中部（z + = 180）桿狀顆粒取向在拉格朗日坐標(biāo)系中的分布（從紅色符號到橙色符號，代表顆粒逐漸從極細(xì)長趨于球形）[38]Fig. 7 Distribution of alignment of particles in the Lagrangian frame near the channel center (z+ = 180)(the symbols from red to orange represent from infinite slender to spherical particles, respectively)[38]

圖8 黏性底層（z +=4）桿狀顆粒取向在拉格朗日坐標(biāo)系中的分布（從紅色符號到橙色符號代表顆粒逐漸從極細(xì)長的桿狀趨于球形）[38]Fig. 8 Distribution of alignment of particles in the Lagrangian frame near the wall (z +=4)( the symbols from red to orange represent from infinite slender to spherical particles, respectively)[38]

圖9 臨界方位角α c與 顆粒形狀δ Λ=1-Λ的關(guān)系（虛線代表斜率1，橫軸從左到右代表顆粒逐漸從極細(xì)長的桿狀趨于球形）[38]Fig. 9 Relationship between critical angle αc andδΛ=1-Λ( black dashed line represents the reference slope δ Λ=1, the values of abscissa from left to right represent from infinite slender to spherical particles, respectively)[38]

4.3 機理分析與模型驗證

4.1節(jié)與4.2節(jié)分別討論了 Λ接近于1的桿狀顆粒在近壁面的區(qū)域依然會與拉格朗日的拉伸方向存在明顯差別。那么產(chǎn)生這種差別的原因是什么？首先，剪切在當(dāng)前問題中扮演著重要的角色，但僅只有剪切這一種因素并不能解釋為什么 當(dāng)Λ→1但Λ≠1的桿狀顆粒在大部分時候會與拉格朗日拉伸方向還具有較好的一致性。因為在只考慮純剪切的情況，根據(jù)Jeffery 方程[26]，不難發(fā)現(xiàn) Λ=1與 Λ≠1兩者顆粒在行為上的迥異差別。對于前者，剪切流中拉格朗日的拉伸方向（ Λ=1）指向流向方向，而壓縮方向會朝著法向。但對于后者，桿狀顆粒與碟狀顆粒均服從Jeffery軌跡在空間內(nèi)周期轉(zhuǎn)動[26]。實際上，對于壁湍流，在壁面附近除了強剪切也存在比較明顯的湍流脈動，湍流脈動引入的噪聲是否會使得 Λ→1但Λ≠1的顆粒在大部分時候還是會與拉格朗日拉伸方向具有較好的一致性？因此，在崔智文等[38]的研究中，速度梯度的脈動被引入。他們認(rèn)為槽道不同區(qū)域的平均剪切與速度梯度的湍流脈動之比(s/D)是影響顆粒取向行為的主要因素，其中s是平均剪切，D表征速度梯度的脈動影響強度。當(dāng)s/D?1時，剪切的作用強于脈動，此時非常小的形狀差異 δΛ也會使得桿狀顆粒與拉格朗日拉伸方向存在明顯的差異；當(dāng)s/D?1時，剪切的作用遠小于脈動，此時比較大的形狀差異 δΛ也會使得桿狀顆粒與拉格朗日拉伸方向基本上表現(xiàn)出較好的一致性。

為了進一步驗證該想法，崔智文等[38]基于Jeffery方程建立了二維簡化的顆粒取向模型，并從方程角度出發(fā)分析在存在強剪切與速度梯度脈動下顆粒轉(zhuǎn)動周期隨形狀的變化，即：

其中，φ是顆粒回轉(zhuǎn)軸在x-z剪切平面內(nèi)的投影與x軸方向的夾角，s是平均剪切， ηφ是隨機脈動。利用文獻[39]中的方法，結(jié)合Fokker-Planck 隨機過程模型，最終可以得到顆粒轉(zhuǎn)動周期統(tǒng)計值的表達式為[38]：

圖10為 顆 粒 平 均 翻 轉(zhuǎn) 時 間 〈τ〉s與 顆 粒 形 狀δΛ=1-Λ的 關(guān)系（ 0＜z+＜10），橫軸從左到右代表顆粒形狀從桿狀逐漸變?yōu)榍蛐?，圖中的符號為直接數(shù)值模擬結(jié)果，虛線為簡化模型的結(jié)果，圖片來自文獻[38]。圖10展示的是理論模型與直接數(shù)值模擬結(jié)果的對比，其中直接數(shù)值模擬結(jié)果通過統(tǒng)計顆粒在所有時刻均停留在范圍 0＜z+＜10的軌跡線片段上的翻滾周期。理論模型中的s由當(dāng)前層的平均剪切得到，而D通過計算拉格朗日軌跡線上的速度梯度自相關(guān)得到。通過圖10可知，理論與數(shù)值實驗結(jié)果的符合度高，進而驗證了他們的想法，即平均剪切與速度梯度的脈動的比值影響顆粒與拉格朗日拉伸方向的相關(guān)性。同時，通過理論模型可知，當(dāng) δΛ→0時，

圖10展示了顆粒的轉(zhuǎn)動周期隨形狀變化 δΛ也存在明顯的平臺區(qū)與冪律區(qū)，而平臺區(qū)的含義代表著顆粒的轉(zhuǎn)動行為基本不隨形狀 δΛ改變，冪律區(qū)則是顆粒周期符合Jeffery理論。這說明當(dāng)顆粒形狀變化進入到冪律區(qū)時，顆粒的行為會產(chǎn)生明顯的變化。同時，區(qū)分平臺區(qū)與冪律區(qū)的臨界值滿足當(dāng)時， δΛc→0， 這意味著非常小的 δΛ也會使得顆粒與Λ=1的情況產(chǎn)生明顯的變化，這與之前分析的結(jié)果一致。

圖10 顆粒平均翻轉(zhuǎn)時間〈 τ〉s與 顆粒形狀δΛ=1-Λ的關(guān)系（0 ＜z+＜10）[38]Fig. 10 Relationship between average period of tumbling of particles and the difference of shape factors of particles δΛ=1-Λ within 0 ＜z+＜10[38]

5 結(jié) 論

本文回顧了非球形顆粒在壁湍流中的取向行為的直接數(shù)值模擬研究。在忽略顆粒慣性的作用下，近壁區(qū)的桿狀顆粒傾向性地朝著流向方向而碟狀顆粒則傾向性地朝著壁面法向方向，且顆粒的傾向性取向行為的程度隨顆粒偏離球形的程度的增加而增強。而非球形顆粒之所以會在壁面具有傾向性取向，是因為偏離球形程度越大的顆粒與流體的拉格朗日拉伸與壓縮方向具有非常強的相關(guān)性。同時，研究工作也發(fā)現(xiàn)即使非球形顆粒的形狀因子 |Λ|→1，但是該類顆粒的取向與拉格朗日拉伸或壓縮方向(| Λ|=1)的在壁面附近還是會存在較大的行為差異。在以桿狀顆粒為例的研究中發(fā)現(xiàn)，桿狀顆粒取向與流體拉格朗日拉伸方向的相關(guān)性隨顆粒形狀的影響規(guī)律在壁面和槽道中心位置是類似的，但不同點在于壁面附近較大的平均剪切與速度梯度脈動之比s/D會使得顆粒取向與流體拉格朗日拉伸方向產(chǎn)生明顯的差異。基于壁面湍流的特點，建立了以平均剪切與速度梯度脈動之比的顆粒轉(zhuǎn)動周期預(yù)測的二維簡化模型。模型的建立與驗證進一步揭示了平均剪切與速度梯度脈動在顆粒取向趨近于拉格朗日拉伸與壓縮方向過程中的重要作用。

在未來的研究工作中，可以考慮非球形顆粒取向行為與湍流中的相干結(jié)構(gòu)之間的相互作用關(guān)系，并將其應(yīng)用在纖維減阻控制的研究之中，進而加深對湍流減阻控制作用的理解。此外，真實流動中非球形顆粒隨著形狀的不規(guī)則性增強，易變形特性會逐漸體現(xiàn)出來，而此時模擬計算中顆粒的柔性需要考慮，所以關(guān)于柔性顆粒在湍流中的行為研究也是該領(lǐng)域可以進一步探索的問題之一。

致謝：本文涉及到工作是作者趙立豪在挪威科技大學(xué)以及回國后在清華大學(xué)的研究內(nèi)容，全部已公開發(fā)表。本工作得到了挪威研究理事會、國家自然科學(xué)基金以及清華大學(xué)國強研究院項目的支持。