單個中性球形顆粒在三維方腔中的運動

李 洋，夏振華

(浙江大學 航空航天學院 流體工程研究所，杭州 310027)

0 引 言

含顆粒多相流現象廣泛存在于自然界和工業生產中，例如流化床、煤燃燒、種子干燥技術、粉體懸浮液輸運、顆粒分離[1-2]、巖土礦石開采中的多孔介質流動等[3]。理解這種類型的多相流對于基礎研究和工程應用都是至關重要的。

很多學者對簡單流動中的顆粒運動進行了大量的研究，諸如顆粒在管道中的沉降[4-5]、顆粒在庫特流中的運動[6-8]和顆粒在管道中的慣性遷移[9-12]。后者稱為慣性聚焦現象， 其廣泛地應用于微通道設備中有限尺寸顆粒的控制和分離[13-14]。在與管道相連的方腔內，Khojah 等[15]發現大（小）顆粒在高雷諾數時向外（內）遷移，而低雷諾數時大（小）顆粒的軌跡靠近（遠離）渦中心。Haddadi 等[16]和Jiang 等[17]也考慮了不同縱橫比的方腔的工況，這些工況也可以用來實現顆粒的分離。

學者們采用實驗方法[18-21]和數值方法[22-27]對頂蓋驅動方腔流進行研究。方腔幾何形狀簡單，但是里面的流動表現出豐富的物理特征，被認為是數值驗證的基準算例。關于無懸浮顆粒的頂蓋驅動方腔流的更多信息，可參閱Shankar 和Deshpande 的綜述[28]。

所有邊界上的壁面約束使得頂蓋驅動方腔流動不同于簡單的槽道流動、管道流動。關于含顆粒的槽道流、管道流的研究已經發表了很多，但是對含顆粒的頂蓋驅動方腔流的研究卻很少。Tsorng 等[29]采用立體成像方法實驗研究了中性懸浮顆粒在頂蓋驅動方腔中的長時間運動軌跡（塑料球形顆粒直徑為3 mm，方腔高度為10 cm）：顆粒在方腔一側做螺旋運動；此外，顆粒的長時間軌跡似乎沿著內部循環模式的優先路徑聚集。后來，Tsorng 等[30]對這項工作進行擴展，研究了剛性顆粒與流體密度比、雷諾數之間的影響關系。隨著雷諾數和斯托克斯數的增加，顆粒軌道向方腔外圍靠近。他們將顆粒的行為解釋為旋轉誘導的力迫使顆粒向優先路徑運動，這與泊肅葉流中的Segre-Silberberg 效應類似[9]。之后，Kosinski 等[1]給出了顆粒團在二維方腔中運動的數值結果，其中雷諾數為1 000，采用雙向耦合方法處理固體點顆粒與流體的相互作用。結果表明顆粒傾向于向壁面遷移，且斯托克斯數越大，遷移行為越強烈。Sidik 和Attarzadeh[2]數值模擬了不同雷諾數下二維方腔內固體點顆粒的運動。他們將模擬結果與之前的實驗結果、數值結果進行比較，以證明三次插值法的應用價值。Safdari 和Kim[31-32]采用格子Boltzmann 方法耦合拉格朗日點顆粒追蹤法，數值研究了含顆粒的三維方腔流動。他們觀察到顆粒的運動軌跡主要取決于顆粒大小、方腔內的渦旋行為和顆粒密度。最近，Hu[33]利用格子Boltzmann 方法詳細研究了單個中性橢球顆粒在二維方腔中的運動。他指出，極限環和顆粒初始朝向與位置無關，而且極限環隨著粒徑的增大向方腔中心收縮，隨著雷諾數的增加向右下角遷移。

綜上所述，有限尺寸顆粒在三維頂蓋驅動流中的動力學尚未得到充分理解。本文采用多松弛格子Boltzmann 方法耦合插值反彈格式模擬了單個中性球形顆粒在三維方腔中的運動。展向的邊界是弱受限的對稱邊界或強受限的固體壁面。著重考慮了初始位置、顆粒大小以及雷諾數對球形顆粒運動的影響。顆粒方腔流的研究能夠加深對顆粒在方管方腔中分離的理解。

1 數學模型

1.1 格子Boltzmann 方法

格子Boltzmann 方法是求解Navier-Stokes 方程的方法之一[34-37]。為了保證數值穩定和參數靈活[38-41]，本文采用了多松弛的格子Boltzmann 模型。不可壓流體的D3Q27 模型對應的顆粒分布函數的演化方程為：

1.2 顆粒的運動

顆粒的平動和轉動由牛頓第二定律和歐拉方程控制：

2 問題設置和數值驗證

2.1 問題設置

球形顆粒在頂蓋驅動方腔中的運動如圖1 所示，一個半徑為R的中性懸浮球形顆粒浸沒在高度為S的立方體方腔中。這里，中性顆粒是指密度等于流體密度的顆粒。方腔的上壁面以恒定的速度U沿著y方向運動，其底壁面和沿y方向的兩個壁面均是靜止的。在z方向采用了兩種不同的邊界條件。對于第一種情況（P0，下稱為展向弱受限情形），展向壁面的限制被弱化，在展向兩側我們使用對稱邊界條件；對于第二種情況（P1，下稱為展向強受限情形），在展向固壁使用無滑移邊界條件。兩個相關的無量綱參數雷諾數Re和顆粒斯托克斯數St，分別定義為：

圖1 球形顆粒在頂蓋驅動方腔中運動的示意圖Fig. 1 Sketch of dynamics of a spherical particle in a lid-driven cubic cavity flow

2.2 數值驗證

為了驗證目前程序的準確性，首先模擬了Re=100、400 和1 000 下的單相頂蓋驅動流，網格為96×96×96。不同雷諾數下，z= 0.5 平面上沿著水平中心線的速度u、沿著豎直中心線的速度v和Ku 等[27]結果的比較分別展示在圖2（a）和圖2（b）中，可以發現模擬結果與Ku 等的結果均吻合良好。

圖2 P1 情況下，z = 0.5 平面上，水平中心線上速度u以及豎直中心線上速度v 的分布Fig. 2 In P1 case, the distribution of u-velocity on the horizontal center line and v-velocity on the vertical center line on the plane z = 0.5

接下來我們對顆粒在方腔中的運動進行了網格無關性測試。在模擬過程中，首先得到穩定的單相流場；然后球形顆粒加入到方腔中特定的位置，并允許其自由運動。對于P0 情況，兩套網格被用來模擬Re= 1 000 下 的 顆 粒 方 腔 流。第 一 組 的 網 格 為96×96×96，其中顆粒半徑為6 個網格，初始網格位置為（77，50，48）；第二組的網格為144×144×144，相應的顆粒半徑為9 個網格，初始位置為（115.5，75，72）。對于顆粒在方腔中的運動，顆粒最終限制在展向z=0.5 平面的極限環上運動。如圖3 所示，兩組模擬下的極限環軌跡很好地重合在一起。此外，我們也發現兩套網格下顆粒受到的力和力矩也能夠較好地吻合。這些結果表明，96×96×96 的網格足以滿足Re= 1 000下的球形顆粒運動的計算精度。在接下來的模擬中，96×96×96 的網格被采用。

圖3 P0 情況下，兩套網格的顆粒最終運動軌跡比較Fig. 3 In P0 case, comparison between the final trajectories of particles between two sets of mesh grids

3 展向弱受限的情形

3.1 初始位置的影響

這一部分我們主要考慮了顆粒初始位置對最終運動的影響。對于P0 情況，雷諾數Re= 1 000 和半徑R= 6 下，對應的St= 0.868。顆粒在展向z= 0.5 平面上的不同初始位置釋放，我們模擬了41 組算例，將顆粒的初始位置展示在圖4 中，虛線為無顆粒時的流函數等值線，1、2、3 為三種最終的穩定運動軌跡。初始位于紅色圓形、藍色菱形或黑色三角形位置的顆粒最終分別在極限環3、極限環2、穩定點1 上運動。根據目前的模擬結果，顆粒的初始位置影響著其最終的運動軌跡，這里存在三種不同的運動模態。當顆粒初始放置在渦中心區域時，顆粒最終遷移到穩定點1，即渦中心，此時幾乎以恒定的角速度繞展向做旋轉運動；當顆粒初始放置在外側的藍色菱形點處，其最終沿著極限環2 做周期性運動；而當顆粒初始位置位于紅色圓點時，球形顆粒最后在極限環3 上做周期性運動。

圖4 P0 情況下，z = 0.5 平面上，顆粒從不同位置釋放的相圖Fig. 4 In P0 case, plane z = 0.5, phase diagram of particles released from different positions

根據圖4 的結果，可以發現極限環3 的軌跡更靠近外圍，它的空間位置不同于極限環2。為了進一步說明顆粒在極限環2 和極限環3 上運動的區別，我們在圖5 中展示了初始網格位置分別位于（77，50，48）和（15，50，48）的兩個算例，在相應極限環3 和極限環2 上運動時的合速度大小以及展向的角速度分量大小隨時間的變化規律（為了方便畫圖，時間軸進行了平移，t0不同）。其中，顆粒速度和角速度分別由上壁面速度U和U/S無量綱化。由圖5 可以發現，在兩個極限環上運動的顆粒，它們的速度和角速度在具體數值上體現出明顯差異。另外，它們都呈現周期性變化，極限環3 和極限環2 對應的無量綱周期分別約為6.9 和7.2。這說明位于極限環3 上的顆粒雖然經歷的周長更長，但是其運動速度更大（圖5（a）），可以用更短的時間完成一個周期的運動。鑒于極限環3 和極限環2 的運動特征高度相似，接下來我們僅選取了其中一種進行詳細分析。

圖5 極限環2 和極限環3 上，顆粒運動的合速度和角速度在一段時間內的變化規律Fig. 5 Time evolutions of the particle’s velocity magnitude and the corresponding angular velocity for the limit cycle 2 and limit cycle 3

為了進一步探索顆粒的運動機理，我們以初始網格位置（77，50，48）為例進行分析。圖6（a）展示了角速度隨時間的演化，角速度由U/S標準化，此時 Ωx和Ωy等 于0， Ωz經過一段時間的不規則運動后周期性地變化。這意味著顆粒經過一段時間的不規則運動后到達極限環，最終在極限環3 上運動。此外顆粒動能E隨時間的演化顯示在圖6（b），這里E=0.5mup2+0.5IΩ2，是平動能和轉動能之和；動能E由 0.5mu2max無量綱化，其中umax代表顆粒最大速度。圖6（b）中虛線框的局部放大圖展示了一個周期上的變化，可以發現，動能E存在最小值，對應圖7的d點。

圖6 顆粒角速度和顆粒動能隨時間的演化Fig. 6 Time evolutions of the angular velocity and kinetic energy for the particle

圖7 為顆粒受力與顆粒速度的夾角，以及顆粒受力在其運動軌跡和其運動軌跡垂直方向上的分量隨時間的演化曲線。圖7（a）給出了極限環上的一個周期，這里a、b、c、d、e、f分別代表顆粒受力與速度垂直時的位置。這里指出當 θ大于90°時，顆粒做減速運動；反之顆粒做加速運動。緊接著我們把顆粒受到的總力F分解成沿著運動軌跡方向的分量Fup和垂直于該運動方向的分量Fuv，力Fup關系到顆粒速度大小的改變，而力Fuv控制著顆粒運動的方向，結果展示在圖7（b）中。可以發現Fup明 顯小于Fuv，而且顆粒在c-d和d-e階段經歷了更大的減速和加速運動。

圖7 顆粒受力F 與速度之間的夾角以及相應顆粒在運動方向上的受力分量Fup 和在垂直于運動方向上的受力分量Fuv 隨時間的變化Fig. 7 The angle between particle force F and velocity and the force component Fup in the direction of motion and the force component Fuv in the direction perpendicular to the direction of motion in function of time

為了更加清晰地觀察顆粒在極限環上的運動規律，幾個典型的位置畫在圖8 中，這里實線代表加速過程，虛線代表減速過程。從圖中可以明顯地看到顆粒交替的加速減速運動，這與圖7 中的結果保持一致。球形顆粒上黑色實線用來表征其轉動特征。由于上壁面的剪切作用，顆粒在方腔中平動的同時伴隨著順時針地轉動。

圖8 z = 0.5 平面上，顆粒在極限環3 上一個周期的運動Fig. 8 Orbiting and rotating of a particle at different positions on the limit cycle 3 in plane z = 0.5

此外我們還測試初始網格位置為（70，50，36）的情形，其他參數保持不變。圖9 為Re= 1 000 時，顆粒位置在x-y平面和展向隨時間的演化曲線。我們發現顆粒交替地向上壁面運動或遠離上壁面。同時一個三維的顆粒軌跡圖展示在圖9（d），有趣的是，我們觀察到球形顆粒只在方腔的一側運動。顆粒首先在橫向螺旋運動的伴隨下向內部遷移，然后逐漸向外圍移動。軌跡呈現兩個明顯的特點：x-y平面上的快速遷移和z方向上相對較慢的橫向運動，這和Tsorng等[29-30]的結果相似。他們實驗觀察到球形顆粒在Re= 860 時被限制在方腔的一側，然而顆粒在高雷諾數Re= 3 200 時在兩側往復運動。

圖9 x-y 平面和展向，顆粒位置隨時間的演化和三維視圖Fig. 9 Evolution of x-y plane and spanwise particle positions over time and 3D views

3.2 雷諾數的影響

這一部分研究雷諾數對固體顆粒運動的影響。初始網格位置選取（74，50，48），顆粒半徑設置為6 個網格。和之前一致，在形成單相穩定流場后球形顆粒加入，并允許其自由運動。對于雷諾數Re= 100、400、1 000，相 應 的 斯 托 克 斯 數 分 別 為St= 0.087、0.347、0.868。對于展向弱受限的情形P0，圖10（a～c）描述了不同雷諾數下的極限環軌跡。其中紅色實心圓點代表顆粒初始位置，實線代表極限環軌跡。背景虛線表示z= 0.5 平面上單相時的流線。極限環上顆粒速度的大小用顏色表示，顆粒速度由相應的上壁面速度U無量綱化。如圖10 所示，隨著雷諾數的增大，極限環逐漸靠近外側，逐漸變大。我們知道單相時，隨著流體慣性的增加，主渦中心向方腔中心移動。極限環軌跡也有往方腔外圍移動的趨勢。在不同雷諾數下，方腔內流場結構存在較大差異，流固相互作用影響了極限環軌跡。

圖10 不同雷諾數下的極限環軌跡Fig. 10 Limit cycles at different Re

此外，我們發現顆粒在極限環軌道上交替的加速和減速運動。當Re= 100 時，顆粒最大速度出現在極限環的頂部，靠近移動的上壁面，最小速度出現在左側。與低雷諾數結果不同的是，Re= 400 和1 000 時的最小速度出現在區域對角線附近的極限環的右上角。如圖7 所示，在d點之前，顆粒的減速周期相對較大且較長。

3.3 顆粒大小的影響

這一部分我們主要考慮了顆粒大小的影響。這里初始網格位置選取（74，50，48），當雷諾數取Re=1 000，顆粒半徑R= 6、9、12，對應的斯托克斯數分別為St= 0.868、1.953、3.472，其他參數保持不變。對于展向弱受限的情形，雷諾數Re= 1 000 時，不同尺寸的顆粒對應的極限環軌跡展示在圖11（a）。其中紅色實心圓點表示初始位置，從內到外分別代表R=6、9、12 時的極限環。背景虛線表示z= 0.5 平面上單相時的流線。極限環上顆粒速度的大小用顏色表示，顆粒速度由相應的上壁面速度U無量綱化。由于初始位置位于平面z= 0.5，顆粒僅在該平面內運動，且在穩定軌道上交替地加速和減速運動。此外，隨著St的增加，極限環向外圍移動。在這里，St描述了方腔流動中懸浮固體顆粒的行為。如果St越大，則表示顆粒受慣性的影響較大，繼續沿著原來的軌跡運動；反之則沿著相應的流體流線運動。如圖11（a）所示，極限環軌跡隨著顆粒大小的變化而變化，大致呈橢圓形。當顆粒運動到極限環的上部時，R= 12 時的St較大，導致與流體流線偏離。由于速度較低，它在底部幾乎沿流線移動。

圖11 不同雷諾數下，極限環大小隨顆粒大小的變化規律Fig. 11 Limit cycles for particles of different sizes at different Re

此外，我們還模擬了Re相對較低的情形。這里采用Re= 100，R= 6、12，對應斯托克斯數分別為St=0.087、0.347。其 他 參 數 與Re= 1 000 的 情 況 相 同。圖11（b）給出了不同顆粒半徑下的極限環，其中紅色實心點為初始位置，外側實線和內側實線分別為R=6 和12 的極限環。與高雷諾數結果不同的是，大的顆粒向靠近渦中心的軌道遷移，而小的顆粒向外極限環遷移，這與Khojah 等[15]的結論定性上一致。根據顆粒或細胞的尺寸特性，通過改變相應的流動慣性或腔體尺寸可以在微流控平臺上實現慣性分離[15-17]。

4 展向強受限的情形

4.1 初始位置的影響

這一部分我們主要考慮了展向在強受限情況下初始位置對顆粒最終運動的影響。對于P1 情況，雷諾數Re= 1 000 和半徑R= 6 的參數下，對應的斯托克斯數St= 0.868。顆粒在展向z= 0.5 平面上的不同初始位置釋放，我們模擬了25 組算例，且顆粒的初始位置如圖12 中空心菱形所示，黑色實線表示最終的穩定軌跡。可以發現，顆粒的初始位置對最終的極限環軌跡沒有影響。這與Hu[33,48]的結果一致，他在文章中指出，二維方腔中顆粒的極限環軌跡與長橢球顆粒的初始朝向和初始位置無關。

圖12 P1 情況下，z = 0.5 平面上顆粒從不同位置釋放下的相圖Fig. 12 Phase diagram for a spherical particle inside lid-driven cavity flow in plane z = 0.5 for cases P1

4.2 雷諾數的影響

這一部分研究雷諾數對固體顆粒運動的影響。初始網格位置選取（74，50，48），顆粒半徑設置為6 個網格。對于雷諾數Re= 100、400、1 000，相應的斯托克斯數St= 0.087、0.347、0.868。對于展向強受限的情形P1，圖13 給出了展向z= 0.5 平面上不同雷諾數下的極限環軌跡，其中紅色實心圓點代表顆粒初始位置。發現雷諾數對極限環的拓撲結構有明顯的影響。隨著雷諾數的增加，除了方腔的左上區域外，極限環軌道有向外遷移的趨勢。Hu[48]在文章中指出，隨著雷諾數的增加，左上角的渦旋發展，球形顆粒在二維方腔中的極限環被推向方腔右下角。然而由于展向強受限（P1），極限環的拓撲結構與P0 情形下的結果有很大不同。

圖13 P1 情況下，Re = 100、400、1 000 的極限環軌跡Fig. 13 Limit cycles at Re = 100、400 and 1 000 for the cases P1

4.3 顆粒大小的影響

最后我們主要考慮了顆粒大小對最終運動軌跡的影響。這里初始網格位置選取（74，50，48），當雷諾數取Re= 1 000，顆粒半徑R= 6、12 對應的斯托克斯數St分別為0.868、3.472，其他參數保持不變。對于展向強受限的情形P1，雷諾數Re= 1 000 時，不同尺寸的顆粒對應的極限環軌跡展示在圖14（a）。其中紅色實心圓點表示初始位置。可以發現小顆粒的極限環在外圍，大顆粒的極限環在內側。這種現象和展向不受限情形P0 下Re= 1 000 的結果恰好相反。

此外，我們還模擬了Re相對較低的情形。這里采 用Re= 100，R= 6、12，對 應St數 分 別 為0.087、0.347。其他參數與Re= 1 000 的情況相同。圖14（b）給出了顆粒半徑R= 6 和12 時的極限環，其中紅色實心點為初始位置，外側實線和內側實線分別為R=6 和12 的極限環，虛線為顆粒到達極限環之前的運動軌跡。與高雷諾數Re= 1 000 的結果一致，大的顆粒向靠近渦中心的軌道遷移，而小的顆粒向外極限環遷移。這與展向弱受限情形P0 下的不同，可能是由于展向壁面的約束導致的。

圖14 不同雷諾數下，極限環大小隨顆粒大小的變化規律Fig. 14 The variation of limit cycle size with particle size at different Re

5 結 論

本文采用多松弛格子Boltzmann 方法模擬了單個中性球形顆粒在三維方腔中的運動。首先通過單相頂蓋驅動流和網格無關性測試驗證了程序的準確性。其次主要考慮了方腔展向弱受限（對稱邊界條件）和強受限（固壁無滑移邊界條件）兩種情況，著重研究了初始位置、顆粒大小以及雷諾數的影響，得到以下結論：

對于展向弱受限的情形P0，發現顆粒的初始位置強烈地影響著最終的運動軌跡。根據展向z=0.5 平面的相圖，可以得到其被劃分為三個區域：外層穩定區、內層穩定區以及渦中心區域。通過對顆粒受力的分解，解釋了其在極限環上運動的機理。此外，還詳細介紹了球形顆粒在極限環上的順時針旋轉運動，呈現交替性地加速和減速運動。也模擬了初始位置位于展向z/S= 0.375 平面的情形，觀察到球形顆粒只在方腔的一側運動，伴隨著在展向z方向上相對較慢的運動。隨著雷諾數的增加，顆粒逐漸向外圍靠近，不斷旋轉達到最優路徑，這里最優路徑是指相應參數下對應的極限環軌跡，運動路徑幾乎遵循相應的流線。對于選定的初始位置，我們觀察到大顆粒的極限環在高雷諾數時向外遷移，而在低雷諾數時更靠近渦中心。

對于展向強受限的情形P1，極限環軌跡與顆粒的初始位置無關。隨著雷諾數的增加，除方腔左上角外，極限環軌跡有向外遷移的趨勢。在Re= 1 000 和100 兩種情況下，隨著顆粒尺寸的增大，極限環都會相應縮小。

通過與參考文獻的對比分析，在一定程度上驗證了本文結果，希望目前的結果對管道和空腔流中顆粒和細胞的過濾和分離提供一些參考。