李 青，陳堅(jiān)強(qiáng)，畢 林，袁先旭，*

(1. 中國空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，綿陽 621000；2. 中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽 621000)

0 引 言

顆粒流體兩相流廣泛存在于化工及航空航天等領(lǐng)域，如流態(tài)化反應(yīng)器、航空發(fā)動(dòng)機(jī)和大推力火箭燃燒室等。在這些實(shí)際問題中，對于彌散相，一般顆粒數(shù)目NP從 O(1 02) 到 O(1 09)，顆粒相對攜帶流體的密度 比從 O (1) 到 O(103)，如化工管道中的顆粒懸浮物或者是河流泥沙懸浮顆粒如氣固流化床、航空航天中的顆粒湍流在某些極端條件下，如高溫顆粒壁湍流，密度比可以達(dá)到O(1 04)，顆粒雷諾數(shù)體的密度，μ 是流體的動(dòng)力黏度，Uf是流體在顆粒重心位置的歐拉速度，Vp是顆粒的拉格朗日速度，a是其中 ρf是流顆粒半徑。從 O(1 0-3) 到 (1 03)，從 O(1 0-3)到 O(1 03)…… 特征參數(shù)覆蓋的范圍相當(dāng)廣泛，問題復(fù)雜，物理內(nèi)容豐富。不同的工業(yè)背景下，我們遇到的問題也不同，因此提煉的科學(xué)問題以及解決的方法也不同。比如化工領(lǐng)域關(guān)心傳熱傳質(zhì)的宏觀量、質(zhì)量和熱量增量，用于指導(dǎo)優(yōu)化過程工程設(shè)計(jì)；航空航天領(lǐng)域則關(guān)心阻力、壓力分布和熱流密度分布，用于優(yōu)化各種飛行器的性能；生物醫(yī)學(xué)領(lǐng)域，比如心腦血管系統(tǒng)則關(guān)心顆粒沉積規(guī)律或者慣性聚集規(guī)律，用來預(yù)測血栓或從血液里過濾癌細(xì)胞，從而更好預(yù)防疾病或者設(shè)計(jì)相關(guān)醫(yī)療器械等。一般地，開發(fā)顆粒流體歐拉-歐拉兩相連續(xù)模型是學(xué)術(shù)界的目標(biāo)，從而為解決工業(yè)界問題提供可靠的技術(shù)手段。歐拉模型把彌散相當(dāng)成連續(xù)介質(zhì)考慮，然后根據(jù)具體問題進(jìn)行各種簡化假設(shè)，建立封閉關(guān)系式，然后用顆粒解析或點(diǎn)顆粒的方法把這些關(guān)系式耦合到歐拉框架下（Li[1], Morris[2]）。為了解決這個(gè)問題，我們需要理解顆粒兩相流在不同幾何邊界下的流體動(dòng)力學(xué)機(jī)制。比如Morris[2]只考慮平行剪切流的情況，駐點(diǎn)流情況就無法適用。駐點(diǎn)流的特點(diǎn)是，遠(yuǎn)場速度大，然后隨著流體質(zhì)點(diǎn)運(yùn)動(dòng)到近壁面，流體質(zhì)點(diǎn)減速，動(dòng)能轉(zhuǎn)化為壓力。因此，（以流體質(zhì)點(diǎn)沿著對稱軸法向運(yùn)動(dòng)為例），沿著法向的靠近壁面，流體質(zhì)點(diǎn)速度不斷減小，壓力不斷增大，當(dāng)流體質(zhì)點(diǎn)運(yùn)動(dòng)到駐點(diǎn)的時(shí)候?qū)M足壁面無滑移和無穿透條件，滯止在駐點(diǎn)，此時(shí)駐點(diǎn)壓力最大，流體質(zhì)點(diǎn)動(dòng)能為零；另一方面，駐點(diǎn)的巨大壓力推動(dòng)流體質(zhì)點(diǎn)遠(yuǎn)離駐點(diǎn)，壓力不斷轉(zhuǎn)化為流體質(zhì)點(diǎn)動(dòng)能，隨著流體質(zhì)點(diǎn)不斷遠(yuǎn)離駐點(diǎn)，流體質(zhì)點(diǎn)速度不斷增加，壓力減小。通過分析得知，駐點(diǎn)流的速度場和壓力場是全場非均勻的。

在工程應(yīng)用中，如頁巖氣開采，通過把帶顆粒的非牛頓流體垂直向下注入地表層中，液體壓力會(huì)使得巖石破裂，而懸浮顆粒會(huì)占據(jù)巖石裂縫的空間，使得巖石無法閉合。 這樣，水壓力就可以進(jìn)一步破裂巖石，從而為開采地下的天然氣提供條件。因此近年來基于工程需求的推動(dòng)，駐點(diǎn)流中的顆粒動(dòng)力學(xué)蓬勃發(fā)展。例如 Vigolo[3]用點(diǎn)顆粒模型耦合直接數(shù)值模擬和物理實(shí)驗(yàn)研究T管流中的大密度比顆粒群運(yùn)動(dòng)（圖1）；Li[4-5]用顆粒解析的直接數(shù)值模擬PR-DNS研究一個(gè)中性顆粒在三維軸對稱駐點(diǎn)流對稱軸上的運(yùn)動(dòng)；Manoorkar[6-7]通過物理實(shí)驗(yàn)和直接數(shù)值模擬，研究T管流中有限尺寸中性顆粒群的輸運(yùn)特征；Rallabandi[8]用理論解研究在低雷諾數(shù)條件下，一個(gè)顆粒在駐點(diǎn)流中的運(yùn)動(dòng)。Vigolo[3]發(fā)現(xiàn)，顆粒反彈系數(shù)和顆粒慣性高度相關(guān)，而且由于三維渦效應(yīng)，顆粒反彈系數(shù)分布是彌散狀的。Manoorkar[6-7]也同樣發(fā)現(xiàn)T管流的三維渦效應(yīng)，還發(fā)現(xiàn)顆粒輸運(yùn)慣性分選的特征。Wang[9]通過物理實(shí)驗(yàn)研究了駐點(diǎn)流的三維渦效應(yīng)，他們通過在遠(yuǎn)離壁面的距離施加不同的擾動(dòng)，得到了一個(gè)類似于卡門渦的脫離渦對，這個(gè)渦對在靠近壁面時(shí)始終不撞擊對稱中心（圖2）。工業(yè)問題中的駐點(diǎn)流非常容易出現(xiàn)三維渦效應(yīng)，抑或是流場本身出現(xiàn)三維渦結(jié)構(gòu)（Vigolo[3]），抑或是來流擾動(dòng)導(dǎo)致（Wang[9]）。鑒于此，Li[4-5]采用三維軸對稱的數(shù)值控制域，考慮中性顆粒只在對稱軸上運(yùn)動(dòng)的情況，移除了駐點(diǎn)流三維不對稱效應(yīng)帶來的復(fù)雜性，得出了相似結(jié)論—顆粒慣性主導(dǎo)駐點(diǎn)流中的顆粒動(dòng)力學(xué)，并且得出了更多的進(jìn)一步細(xì)節(jié)。

圖1 T型管流運(yùn)動(dòng)的示意圖Fig. 1 Particles clustering in a T junction

除了Vigolo[3]的大密度比顆粒群課題，大多數(shù)研究的共同點(diǎn)都是：顆粒慣性與流體慣性相當(dāng)St～O(1)，密度比背景流動(dòng)是全場非均勻駐點(diǎn)流。這種情況下，一方面，顆粒慣性和環(huán)境流體慣性相當(dāng)；另一方面，顆粒的有限尺度往往和流體微團(tuán)最小水動(dòng)力尺度相當(dāng)，展現(xiàn)出高度的局部各向異性特征，換句話說，顆粒相與流體進(jìn)行動(dòng)量交換的方式是空間非均勻的[10]（圖3）。因此，常規(guī)的點(diǎn)顆粒模型數(shù)值方法將不再適用，必須采用考慮有限尺寸效應(yīng)的顆粒解析的數(shù)值算法：如有限元、格子玻爾茲曼法LBM[11]， 貼體網(wǎng)格下的有限差分或有限體積[12]，虛擬點(diǎn)法[13]，浸沒邊界方法（Immersed Boundary Method，IBM）[14]。考慮顆粒的有限體積效應(yīng)一直是顆粒流體兩相流的挑戰(zhàn)之一，尤其是在顆粒湍流中，有限體積的顆粒注入各向異性的動(dòng)量擾動(dòng)，會(huì)影響湍流的脈動(dòng)統(tǒng)計(jì)量[15]。因此，為了保證數(shù)值實(shí)驗(yàn)的正確性，Li[4-5]采用了浸沒式邊界方法IBM耦合有限體積流體求解器的數(shù)值算法對駐點(diǎn)流的顆粒動(dòng)力學(xué)進(jìn)行分析。

圖2 駐點(diǎn)流來流擾動(dòng)渦對在近壁面的非對稱特性[9]Fig. 2 The non-axisymmetric characteristic of a perturbed vortex pair in a stagnation point flow[9]

圖3 有限尺寸懸浮顆粒尺度與流體微團(tuán)最小運(yùn)動(dòng)尺度量級相當(dāng)情況下的示意圖Fig. 3 A schematic of a suspended particle with the same order of characteristic length scale as that of the flow

接下來，簡要介紹和回顧Li[4-5]的核心工作。Li特別關(guān)注駐點(diǎn)流中的顆粒動(dòng)力學(xué)機(jī)制問題，這個(gè)領(lǐng)域還罕有研究，但是卻很重要。在實(shí)際工業(yè)問題中，顆粒懸浮的裝置往往是復(fù)雜幾何外形，包括駐點(diǎn)流情況。開發(fā)顆粒兩相流的歐拉-歐拉連續(xù)模型則首先需要對復(fù)雜幾何外形下的顆粒懸浮機(jī)理有清晰的認(rèn)識(shí)。駐點(diǎn)流顆粒懸浮問題是一類特殊和典型的問題，比如常見的化工T型管道。

1 顆粒解析的直接數(shù)值模擬方法：浸沒式邊界方法IBM

方程(1)和方程(2)概括描述了顆粒解析的直接數(shù)值模擬方法IBM。式中 Ω 是顆粒體積，u是流場速度矢量，Vp是顆粒速度矢量，是滿足流體與剛性顆粒邊界無滑移條件的動(dòng)量交換矢量,α是光滑函數(shù)以確保數(shù)值穩(wěn)定性，是在顆粒表面拉格朗日點(diǎn)上流體相和顆粒相的滑移速度矢量，dt是流體求解器的時(shí)間步長[14]。一般來說，IBM方法可以計(jì)算顆?；评字Z數(shù)Rep～O(100)的問題，IBM方法也可以用在復(fù)雜外形計(jì)算湍流問題[16]。

在湍流問題中，隨著Re增加，對數(shù)值算法的網(wǎng)格解析精度要求增加，因?yàn)槎喑叨鹊耐牧鞯淖钚∶}動(dòng)尺度K41隨之減小[17]。值得注意的是，流場雷諾數(shù)和顆?；评字Z數(shù)是兩個(gè)不同的概念。當(dāng)流場雷諾數(shù)Re很高的時(shí)候，顆?；评字Z數(shù)Rep可以很小。一般來說，工業(yè)問題中的高雷諾數(shù)顆粒湍流，Rep～O(100)。所以，IBM 方法可以用來研究高雷諾數(shù)顆粒湍流。

2 單個(gè)中性顆粒在駐點(diǎn)流中的動(dòng)力學(xué)

駐點(diǎn)流的理論模型Hiemenz 流動(dòng)可以在三維或者三維軸對稱的條件下，通過自相似解求得[18]。Hiemenz 流動(dòng)是全場非均勻的，遠(yuǎn)場速度大，近壁面速度小，壓力則是遠(yuǎn)場小，近壁面大，這是因?yàn)榱黧w動(dòng)能轉(zhuǎn)化為壓力的緣故。遠(yuǎn)離壁面，壁面法向距離大于大約三個(gè)Hiemenz 特征長度δ，也叫黏性邊界層厚度，Hiemenz理論解求出的渦量幾乎為零[18]，流場是無黏的，是無黏擠壓流動(dòng)，擠壓率可表示為B，根據(jù)具體問題不同，擠壓率可以橫跨好幾個(gè)數(shù)量級B=O(1)～O(104)。 在化工流動(dòng)中，一般擠壓率B～O(1)，以T 型管流為例，B=U/D，U是垂直管流的平均速度，D是等直徑T型管的直徑。而在航空航天飛行器的頭部駐點(diǎn)流中，擠壓率可以達(dá)到B～O(104)，B=U/D，U是脫體激波后的流體速度，D是脫體激波的距離?？拷诿?，由于壁面的無滑移滯止作用，黏性邊界層會(huì)發(fā)展起來，Hiemenz理論解求出的渦量幾乎為零[18]，Hiemenz流動(dòng)中的黏性邊界層厚度是常數(shù)v是流體的動(dòng)力黏度，Bδ是流體的特征速度。因?yàn)镠iemenz流動(dòng)是全場非均勻的，所以它沒有特征流動(dòng)雷諾數(shù)，只有局部雷諾數(shù)Rer=(r/δ)2，Rez=(z/δ)2，所以，通過局部最大雷諾數(shù)來控制數(shù)值計(jì)算域，防止數(shù)值發(fā)散[4]。防止數(shù)值發(fā)散至少有好幾個(gè)方法，目前已知的一個(gè)是局部最大雷諾數(shù)不能太大，一個(gè)是遠(yuǎn)場采用非均勻粗網(wǎng)格。

圖4表明，Li[4-5]的數(shù)值方法能很好地重構(gòu)一個(gè)Hiemenz流場并與理論解吻合。 圖4 的橫軸分別表示歸一水平速度和法向速度，豎軸表示歸一化的壁面距離，黑色虛線表示Hiemenz流動(dòng)的理論解，從左到右不同彩色實(shí)線表示不同徑向位置切面r/δ=1.56，3，8的DNS計(jì)算的速度剖面結(jié)果。Li[4]在此基礎(chǔ)上研究單個(gè)有限尺寸中性懸浮顆粒在駐點(diǎn)流中的動(dòng)力學(xué)（圖5）。

2.1 與壁面發(fā)生彈性碰撞前

圖6 顯示了不同有限尺寸顆粒受到的無黏環(huán)境壓力的理論解，Ba2是駐點(diǎn)流的特征水動(dòng)力[19]，見式(3～5)。圖6中橫軸表示不同有限尺寸中性懸浮顆粒在駐點(diǎn)流中受到的歸一化無黏環(huán)境壓力的理論解，豎軸表示歸一化的壁面距離，從左到右不同的彩色線表示從小到大的有限尺寸中性懸浮顆粒a/δ=0.8，1.6，2.4，3.2?？梢钥吹剑绻邢蕹叽鐟T性顆粒只受到無黏環(huán)境壓力作用，那么顆粒的運(yùn)動(dòng)導(dǎo)數(shù)和流體微團(tuán)的隨體導(dǎo)數(shù)是等價(jià)的，那這種情況可以認(rèn)為顆粒有示蹤粒子行為[20]。

圖4 DNS流場歸一化速度與駐點(diǎn)流理論解驗(yàn)證Fig. 4 Comparisons of dimensionless velocities between DNS and the theoretical solution

圖5 一個(gè)有限尺寸中性懸浮顆粒在三維軸對稱駐點(diǎn)流中的動(dòng)力學(xué)的顆粒解析的直接數(shù)值模擬PR-DNS示意圖[5]（ω/B表示歸一化標(biāo)量渦量）Fig. 5 A finite-size neutrally-buoyant particle in an axisymmetric stagnation point flow obtained by the particle-resolved DNS[5](ω/B represents thedimensionless vorticity intensity)

圖6 不同有限尺寸顆粒受到的無黏環(huán)境壓力的理論解[4]Fig. 6 The dimensionless theoretical ambient pressure force excerted on different finite size neutrally-buoyant particle[4]

圖7（a）顯示了不同有限尺寸顆粒沿法向?qū)ΨQ軸受到的總水動(dòng)力。遠(yuǎn)離壁面，總水動(dòng)力線性遞減，只受到無黏環(huán)境壓力，呈示蹤粒子行為。靠近壁面，由于壁面黏性邊界層與顆粒相互作用，顆粒受到的總水動(dòng)力不再是無黏環(huán)境壓力，示蹤粒子行為不再成立，見圖7（b）。

圖7 在不同壁面距離條件下，不同慣性顆粒受到的總水動(dòng)力以及速度與流體微團(tuán)速度的比較[4]Fig. 7 The hydrodynamic force on inertial particles at different gaps between the bottom of particles and the wall , and the velocities of inertial particles as a function of the distances between the center of particles and the wall [4]

圖7（a）的橫軸表示不同有限尺寸顆粒在駐點(diǎn)流中受到的歸一化總水動(dòng)力，豎軸表示歸一化的壁面距離，從左到右不同的彩色線表示從小到大的有限尺寸中性懸浮顆粒a/δ=0.8，1.6，2.4，3.2。圖7（b）的橫軸表示不同有限尺寸顆粒在駐點(diǎn)流中的歸一化速度，豎軸表示歸一化的壁面距離，藍(lán)色實(shí)線表示流體示蹤粒子，其他的彩色線與圖7（a）意義相同。 圖7上圖表明靠近壁面，不同慣性顆粒受力趨勢分為兩種情況，一種情況顆粒受力趨于零，另一種趨于發(fā)散。這是因?yàn)榻诿娴念w粒受力主要由潤滑力控制，，當(dāng)顆粒慣性小的時(shí)候，潤滑力使得慣性顆粒動(dòng)能耗散為零V～0，潤滑力也收斂為零Flub～0；另一種情況，潤滑力呈現(xiàn)發(fā)散趨勢，因?yàn)閼T性作用還不能在有限壁面距離內(nèi)被流體的黏性作用完全耗散，顆粒以有限大小速度撞擊壁面，而 ε～0，所以呈發(fā)散趨勢，這里ε 是顆粒與壁面的無量綱壁面距離。 這種情況下，彈性碰撞將會(huì)發(fā)生，純粹描述流體動(dòng)力學(xué)的N-S方程將不再能描述這個(gè)復(fù)雜物理問題，物理問題將變?yōu)榱鞴恬詈蠁栴}，必須考慮固體彈性碰撞的過程。

2.2 與壁面發(fā)生彈性碰撞后

當(dāng)彈性碰撞發(fā)生后，方程（6）描述了顆粒動(dòng)力學(xué)，方程（1）描述了有顆粒擾動(dòng)的流體力學(xué)場。理論上是，求解一個(gè)考慮顆粒相動(dòng)量注入的流體力學(xué)偏微分方程，同時(shí)求解一個(gè)顆粒動(dòng)力學(xué)方程.之前的難點(diǎn)在于如何給出Fc的表達(dá)式。因?yàn)轭w粒彈性碰撞力是物性材料和碰撞時(shí)間的函數(shù)，這里簡要表示為Fc=f(Tc)，Tc是理論碰撞力表達(dá)式的理論碰撞時(shí)間。一般地，當(dāng)顆粒慣性大的時(shí)候 St～O(1000)，顆粒反彈系數(shù)和流場特性無關(guān)，近似等效于顆粒本身的干碰撞系數(shù)[21-23]。但是當(dāng)顆粒慣性小的時(shí)候 St～O(10)，反彈系數(shù)就和顆粒慣性相關(guān)，不再直接等于干碰撞系數(shù)，這是因?yàn)樗畡?dòng)力效應(yīng)黏性耗散了部分顆粒動(dòng)能[21-22]。既然如此，當(dāng)顆粒慣性小的時(shí)候，顆粒碰撞時(shí)間也沒有理由是常數(shù)[24]。Li[5]總結(jié)了前人的實(shí)驗(yàn)和理論結(jié)果（ Zenit[23]和Birwa[24]）后，采用彈性力學(xué)解析解（Goldsmith[25]）來表達(dá)數(shù)值模擬中的Tc，獲得了很好的結(jié)果。

圖8表明，Li[5]的數(shù)值方法能很好地重構(gòu)濕潤碰撞的物理場景，準(zhǔn)確預(yù)測了濕潤碰撞的反彈系數(shù)，并且反應(yīng)濕潤碰撞的不同慣性顆粒的不同碰撞時(shí)間，濕潤碰撞時(shí)間不取決于任意的數(shù)值光滑參數(shù)，數(shù)值上具有魯棒性。 圖8（a）的橫軸表示歸一化的不同顆粒慣性，豎軸表示歸一化的反彈速度，黑白符號代表實(shí)驗(yàn)數(shù)據(jù)（Joseph[21]和 Gondret[22]），彩色符號代表不同數(shù)值參數(shù)下的直接數(shù)值模擬DNS結(jié)果。圖8（b）的橫軸表示歸一化的不同顆粒慣性，豎軸表示歸一化的顆粒濕潤碰撞時(shí)間，插圖的豎軸表示歸一化的干碰撞時(shí)間，彩色符號代表不同數(shù)值參數(shù)下的直接數(shù)值模擬DNS結(jié)果。Nc表示數(shù)值光滑參數(shù)，η表示無量綱歸一化粗糙度，Stc表 示不同粗糙度 η條件下的顆粒反彈發(fā)生的臨界顆粒慣性。在化工領(lǐng)域，顆粒碰撞的典型慣性一般為是顆粒在靜止環(huán)境流體中的定常沉降速度。所以Li[5]的數(shù)值方法是有針對性的。

圖8 PR-DNS數(shù)值模擬的重力沉降問題的反彈系數(shù)和濕碰撞時(shí)間[5]Fig. 8 PR-DNS results of rebound coefficients and wet collision time for settling problem[5]

Li[5]給出了不同慣性的有限尺寸中性顆粒的反彈系數(shù)，數(shù)值實(shí)驗(yàn)表明，中性顆粒在駐點(diǎn)流中的動(dòng)力學(xué)與沉降顆粒類似，都是依賴于顆粒慣性；數(shù)值方法上，濕潤碰撞時(shí)間反應(yīng)濕潤碰撞的不同慣性顆粒的不同碰撞時(shí)間，濕潤碰撞時(shí)間不取決于任意的數(shù)值光滑參數(shù)，數(shù)值上具有魯棒性，見圖9。這里Nc表示數(shù)值光滑參數(shù)，η表示無量綱歸一化粗糙度，ac表示不同粗糙度 η條件下的顆粒反彈發(fā)生的臨界顆粒半徑。圖9（a）的橫軸表示歸一化的不同顆粒慣性，豎軸表示歸一化的反彈速度，彩色符號代表不同數(shù)值參數(shù)下的直接數(shù)值模擬DNS結(jié)果。圖9（b）的橫軸表示歸一化的不同顆粒慣性，豎軸表示歸一化的顆粒濕潤碰撞時(shí)間，彩色符號代表不同數(shù)值參數(shù)下的直接數(shù)值模擬結(jié)果。 我們發(fā)現(xiàn)，中性顆粒反彈的閾值、轉(zhuǎn)捩點(diǎn)在a/δ～1處，這是合乎常理的。因?yàn)?，一開始Li[4]把數(shù)值實(shí)驗(yàn)設(shè)置為中性顆粒，移除了兩相密度比對顆粒慣性的影響，換句話說，顆粒慣性完全來自于顆粒本身的有限尺寸。如果顆粒有限尺寸與邊界層厚度相當(dāng)，那么顆粒慣性動(dòng)能將會(huì)完全被黏性邊界層耗散掉，見圖7、圖9和圖10。如果顆粒有限尺寸比邊界層厚度大很多，那么顆粒慣性將無法在有限截?cái)喑叨鹊谋诿婢嚯x內(nèi)被黏性邊界層完全耗散（具體是潤滑力主導(dǎo)這個(gè)過程），那么慣性顆粒將會(huì)以有限速度撞擊壁面，見圖7，彈性碰撞和反彈將會(huì)發(fā)生，見圖9?？傊?，單個(gè)有限尺寸中性顆粒在駐點(diǎn)流中的動(dòng)力學(xué)將完全由顆粒慣性a/δ控制[4]。

圖10 不同有限尺寸中性顆粒與黏性邊界層厚度對比示意圖，a）和b）分別表示不同的有限尺寸中性懸浮顆粒在三維軸對稱駐點(diǎn)流近壁面對流場的影響[4]Fig. 10 A schematic of surrounding flow fields around neutrallybuoyant particles with different sizes. Dashed lines represent the boundary layer edge [4]

3 一對中性顆粒在駐點(diǎn)流中的動(dòng)力學(xué)

實(shí)際問題中，顆粒往往是一群的，它們聚集或者分散，因此研究一對中性顆粒的行為將為我們研究大量彌散顆粒在駐點(diǎn)流中的動(dòng)力學(xué)打下研究基礎(chǔ)（見圖11（a））。在確保網(wǎng)格無關(guān)和時(shí)間步無關(guān)的前提下，我們測試了不同的慣性顆粒和不同的初始顆粒間距，確保圖11（b）的結(jié)果是可重復(fù)的[5]。 圖11（b）的橫軸表示歸一化的時(shí)間，豎軸表示不同顆粒的歸一化壁面距離，藍(lán)色虛線代表單個(gè)中性顆粒，藍(lán)色實(shí)線代表兩倍于藍(lán)色虛線的顆粒半徑的中性顆粒，紅黑線代表分別代表一對顆粒中的低位和高位顆粒，紅色實(shí)線代表與藍(lán)色虛線半徑相同的低位顆粒，黑色實(shí)線代表與藍(lán)色虛線半徑相同的高位顆粒，綠色符號表示低位顆粒的無壁面接觸的反彈點(diǎn)。 我們發(fā)現(xiàn)，一對顆粒中的低位顆粒（紅色軌跡）可以不接觸壁面發(fā)生反彈，見圖11（b）和插圖里的綠色符號，而單個(gè)顆粒卻沒有這種行為（藍(lán)色虛線軌跡）。這很反常，僅憑借水動(dòng)力效應(yīng)就能把一個(gè)慣性顆粒彈起嗎？

圖11 駐點(diǎn)流中的顆粒非接觸碰撞反彈現(xiàn)象[5]Fig. 11 Contactless bouncing phenomenon of particle in a stagnation point flow[5]

圖12解釋了這種反常的反彈。圖12的水平軸表示歸一化時(shí)間。圖12（a）的豎軸表示不同顆粒的歸一化壁面距離，圖12（b）的豎軸表示不同顆粒受到的總水動(dòng)力，圖12（c）的豎軸表示不同顆粒受到的無黏環(huán)境壓力，圖12（d）的豎軸表示不同顆粒受到的總水動(dòng)力減去無黏環(huán)境壓力的凈效應(yīng)。藍(lán)色虛線代表單個(gè)中性顆粒，紅黑線代表分別代表一對顆粒中的低位和高位顆粒，紅色實(shí)線代表與藍(lán)色虛線半徑相同的低位顆粒，黑色實(shí)線代表與藍(lán)色虛線半徑相同的高位顆粒。 簡單來說，就是一對顆粒的情況下，高位顆粒替低位顆粒阻擋了來自上方的駐點(diǎn)流向下擠壓效應(yīng)[26],而低位顆粒沒有物體替它阻擋。駐點(diǎn)流的非均勻性會(huì)產(chǎn)生一個(gè)向上托舉的無黏環(huán)境壓力[19]，見方程（3），見圖6。圖12也通過顆粒解析的直接數(shù)值模擬的結(jié)果表明，低位顆粒的總水動(dòng)力是向上的，而單個(gè)顆粒則是向下的。無黏環(huán)境壓力始終存在于駐點(diǎn)流，由于一對顆粒的高位顆粒的庇護(hù)作用，使得低位顆?？梢员粺o黏環(huán)境壓力托舉起來[5]。

Li[26]還發(fā)現(xiàn)慣性小的顆粒會(huì)抱團(tuán)撞擊壁面，慣性大的則無法抱團(tuán)，見圖13，水平軸表示歸一化的顆粒慣性；左，右豎軸分別表示單個(gè)顆粒和一對顆粒中的低位顆粒的歸一化反彈系數(shù)。Li[26]認(rèn)為這是駐點(diǎn)流的擠壓效應(yīng)與顆粒慣性競爭的結(jié)果，當(dāng)顆粒慣性小的時(shí)候，擠壓效應(yīng)大于顆粒慣性，顆粒間相互碰撞后無法彈開，被擠壓流牢牢抱緊在一起；慣性大則不會(huì)。Li[26]歸一化了抱團(tuán)顆粒對的慣性后發(fā)現(xiàn)，抱團(tuán)顆粒對的歸一化反彈系數(shù)曲線和單個(gè)顆粒的反彈系數(shù)曲線是重合。這再一次印證了Li[4]的結(jié)論，駐點(diǎn)流中的顆粒動(dòng)力學(xué)由顆粒慣性控制。

圖12 單個(gè)和一對有限尺寸中性顆粒在駐點(diǎn)流中的受力以及各力分量比較示意圖[5]Fig. 12 A illustration of the dynamics of a single and a pair of neutrally-buoyant particles in a stagnation point flow[5]

圖13 單個(gè)（藍(lán)色）和一對有限尺寸中性顆粒（紅色）在駐點(diǎn)流中的反彈系數(shù)比較示意圖[26]Fig. 13 The dimensionless rebound velocities of a single (blue)and a pair of (red) isolated neutrally-buoyant particles[26]

4 結(jié)論與展望

駐點(diǎn)流中的顆粒動(dòng)力學(xué)與非駐點(diǎn)流顆粒動(dòng)力學(xué)非常不同。研究通過選取有限尺寸中性懸浮顆粒在三維軸對稱駐點(diǎn)流對稱軸上運(yùn)動(dòng)這個(gè)基本模型，移除了真實(shí)來流擾動(dòng)或者數(shù)值擾動(dòng)的影響，移除了密度比對顆粒慣性的影響，這樣顆粒慣性完全來自于顆粒有限尺寸，顆粒動(dòng)力學(xué)將局限于法向，顆粒與環(huán)境流體和壁面的相互作用，不需要考慮來流擾動(dòng)導(dǎo)致的非定常效應(yīng)，不需要考慮顆粒不在對稱軸上導(dǎo)致的剪切誘導(dǎo)升力。

遠(yuǎn)離壁面，中性顆粒會(huì)像示蹤流體微團(tuán)一樣運(yùn)動(dòng)。通過數(shù)學(xué)推導(dǎo)，研究得出，顆粒動(dòng)力方程的牛頓慣性項(xiàng)等于流體微團(tuán)的隨體導(dǎo)數(shù)。這在數(shù)學(xué)上解釋了：為什么有限尺寸慣性顆粒與數(shù)學(xué)上無窮小的流體微團(tuán)示蹤粒子的運(yùn)動(dòng)軌跡重合。

靠近壁面，顆粒慣性的不同將導(dǎo)致顆粒偏移流體示蹤顆粒軌跡的幅值不同: 慣性大的顆粒，將會(huì)誘導(dǎo)出大的潤滑力，反之亦然。但是潤滑力是數(shù)學(xué)理想模型。當(dāng)顆粒慣性小，黏性潤滑力足以耗散掉所有的顆粒動(dòng)能，顆粒將會(huì)減速至零并停留在壁面；當(dāng)顆粒慣性大，在撞擊有限尺度的粗糙度壁面前，潤滑力不足以耗散掉顆粒動(dòng)能，顆粒將會(huì)以有限速度撞擊壁面，彈性碰撞將會(huì)發(fā)生，這個(gè)有限速度與顆粒慣性成正比。當(dāng)彈性碰撞發(fā)生后，流體力學(xué)問題變成了流固耦合問題。因此我們耦合一個(gè)固體彈性力學(xué)的解析表達(dá)式到碰撞模型里，與流體求解器耦合求解，與經(jīng)典沉降實(shí)驗(yàn)數(shù)據(jù)吻合。

研究發(fā)現(xiàn)，單個(gè)中性顆粒在駐點(diǎn)流中的反彈系數(shù)與沉降顆粒類似，與顆粒慣性相關(guān)。但是當(dāng)研究到一對中性顆粒的時(shí)候，反常的現(xiàn)象發(fā)生了。首先，根據(jù)顆粒慣性不同，顆粒對有可能會(huì)抱團(tuán)，這是駐點(diǎn)流擠壓效應(yīng)與顆粒慣性相互競爭的結(jié)果；其次，當(dāng)顆粒對以抱團(tuán)撞擊壁面的時(shí)候，顆粒的反彈系數(shù)與單個(gè)顆粒反彈系數(shù)是重合的；最后，當(dāng)顆粒對不抱團(tuán)的時(shí)候，低位顆粒可以不撞擊壁面發(fā)生反彈。

研究闡明了駐點(diǎn)流中的有限尺寸中性顆粒動(dòng)力學(xué)完全由顆粒慣性/有限尺寸控制的機(jī)理,為進(jìn)一步研究駐點(diǎn)流中的顆粒動(dòng)力學(xué)鋪墊了基礎(chǔ)。對于某些工業(yè)應(yīng)用中的駐點(diǎn)流中的顆粒動(dòng)力學(xué)問題，顆粒運(yùn)動(dòng)軌跡往往不是在對稱軸上的，那么隨之而來的平行加速流帶來的剪切誘導(dǎo)升力就必須要考慮；此外，顆粒與攜帶流體的密度比也不總是1，顆粒的慣性將不再由有限尺寸單獨(dú)控制，顆粒相對流體的密度比效應(yīng)也需要考慮。未來我們將針對這些真實(shí)問題帶來的復(fù)雜效應(yīng)，進(jìn)行更深入的基礎(chǔ)研究。

致謝感謝國家重點(diǎn)研發(fā)計(jì)劃（2019YFA0405200）經(jīng)費(fèi)的支持，感謝圖盧茲大學(xué)Abbas 教授、Climent教授、Magnaudet教授和紐約城市大學(xué)Morris教授在作者博士工作期間給予的幫助和支持。

APPENDIX: NOMENTCLATURE

ρfcarried fluid density, kg/m3

ρpparticle density, kg/m3

Npparticle number, 1

aparticle radius, m

dparticle diameter, m

Ufunperturbed fluid velocity at particle center, m/s

Vpparticle velocity vector at particle center, m/s

Vnormal particle velocity at particle center, m/s

gnormal gravity constant, m2/s

qparticle flux over unit surface, m2/s

μdynamics viscosity,Pa·s

v kinematic viscosity, m2/s

Bflow strain rate, 1/s

δ Hiemenz flow character length scale, m

Stparticle Stokes number, represents surrounding particle inertia effect in comparing viscous fluid effect, 1