劉 依， 聶德明

（中國計(jì)量大學(xué) 計(jì)量測(cè)試工程學(xué)院，杭州310018）

1 引 言

顆粒的沉降不僅普遍存在于自然界，而且具有廣泛的工業(yè)生產(chǎn)背景，是流體力學(xué)的一個(gè)經(jīng)典問題。影響沉降的因素很多，如顆粒粒徑、流體粘度、顆粒與流體的密度差、顆粒形狀和顆粒濃度等，因此對(duì)顆粒沉降的研究具有一定的復(fù)雜性。從計(jì)算的角度來看，目前主要通過直接數(shù)值模擬方法對(duì)顆粒沉降現(xiàn)象及機(jī)理進(jìn)行研究。

當(dāng)多個(gè)顆粒在流體中一起運(yùn)動(dòng)時(shí)，其相互作用無論對(duì)于顆粒還是流體而言，都至關(guān)重要。尤其在流體的慣性不能忽略的情況下，顆粒之間的相互作用往往會(huì)對(duì)顆粒的運(yùn)動(dòng)形態(tài)產(chǎn)生顯著的影響。DKT（drafting－kissing－tumbling）就是最典型的例子。該現(xiàn)象最早由Fortes等［1］通過實(shí)驗(yàn)觀察到。Feng等［2］則采用有限元方法模擬了兩個(gè)圓形顆粒的沉降，其顆粒雷諾數(shù)設(shè)定為70，其結(jié)果重現(xiàn)了DKT現(xiàn)象。此后，Aidun等［3］通過LBM研究了雙顆粒在無限長通道中的沉降特性，其結(jié)果更為復(fù)雜而豐富，研究發(fā)現(xiàn)，在較低雷諾數(shù)下顆粒的沉降呈現(xiàn)周期性的振動(dòng)，隨著雷諾數(shù)的增大，顆粒的運(yùn)動(dòng)會(huì)出現(xiàn)分岔現(xiàn)象，而進(jìn)一步增大雷諾數(shù)則會(huì)轉(zhuǎn)變?yōu)殡S機(jī)狀態(tài)。Verjus等［4］對(duì) Aidun等［3］的結(jié)果進(jìn)行了補(bǔ)充。文獻(xiàn)［5－7］研究了中等雷諾數(shù)下顆粒的相互作用對(duì)沉降的影響。Yacoubi等［5］發(fā)現(xiàn)顆粒個(gè)數(shù)的奇偶性與沉降形態(tài)的凹凸形狀息息相關(guān)，而Nie等［6，7］則報(bào)道了多個(gè)顆粒在沉降中發(fā)生 DKT的選擇性組合以及分組現(xiàn)象。Wang等［8］采用LBM對(duì)大小不同的顆粒沉降進(jìn)行了細(xì)致研究，重點(diǎn)關(guān)注顆粒之間的尺寸差異對(duì)DKT的影響。邵雪明等［9］采用拉格朗日乘子／虛擬區(qū)域方法，模擬了尺寸不同的雙顆粒的沉降過程，結(jié)果表明，只有當(dāng)顆粒之間尺度差異很小時(shí)，才會(huì)重復(fù)進(jìn)行DKT。曾卓雄等［10］采用數(shù)值方法研究了湍流中顆粒碰撞與湍流脈動(dòng)的關(guān)系。胡平等［11］采用LBM研究了傾斜通道內(nèi)的多顆粒沉降特性，發(fā)現(xiàn)當(dāng)傾斜角處于特定的范圍內(nèi)時(shí)，會(huì)出現(xiàn)后一個(gè)顆粒超越前一個(gè)顆粒的現(xiàn)象。陳榮前等［12］采用LBM對(duì)雙顆粒沉降的動(dòng)力學(xué)進(jìn)行了分析，揭示了中等雷諾數(shù)下顆粒沉降的橫向受力特性。

綜上所述，顆粒之間的相互作用對(duì)顆粒的運(yùn)動(dòng)起關(guān)鍵作用。Aidun等［3］發(fā)現(xiàn)即使對(duì)于雙顆粒沉降的簡(jiǎn)單系統(tǒng)，其呈現(xiàn)的動(dòng)力學(xué)特性也極其復(fù)雜而豐富。在沉降的諸多因素中，顆粒的密度決定了其慣性，當(dāng)密度不同的顆粒置于同一個(gè)流場(chǎng)中時(shí)，其相互作用必定會(huì)受到密度差異的影響，可以預(yù)測(cè)，這種影響在流體慣性存在下會(huì)變得極為顯著，而輕重顆粒的沉降形態(tài)也必定異于相同顆粒的結(jié)果。這種差異以何種方式體現(xiàn)，重顆粒是否一定會(huì)比輕顆粒沉降快，什么條件下重顆粒才會(huì)離開輕顆粒，這在以往的研究中鮮有報(bào)道，因此要回答這些問題還必須通過更為細(xì)致而具體的研究，這是本文工作的出發(fā)點(diǎn)。

2 模型

2．1 物理模型

兩個(gè)顆粒并列放置于無限長的豎直通道中，通道寬度為L，顆粒直徑為d，顆粒密度分別為ρp（左）和（右）。引入一個(gè)無量綱參數(shù)k來描述顆粒之間密度的差異，即k＝（－ρp）／ρf，其中ρf為流體密度，計(jì)算中固定為ρp＝1．5與ρf＝1，且k＞0。本文雷諾數(shù)定義為Re＝U0d／ν，其中ν為流體的運(yùn)動(dòng)粘性系數(shù)，U0為特征速度，可表示為

式中g(shù)為重力加速度，速度選取的規(guī)則與文獻(xiàn)［5－7］相同，即假定顆粒沉降所受阻力與重力和浮力平衡，且阻力系數(shù)為1。其他參數(shù)設(shè)置如下，d＝30，L＝5d，顆粒中心間距固定為2d。在物理單位中，d＝0．15cm，ρp＝1．5g／cm3，ρf＝1g／cm3。網(wǎng)格間距和時(shí)間步長分別為Δx＝5×10－6m和Δt＝5×10－6s。為了模擬無限長通道的邊界，計(jì)算中重顆粒始終距上游邊界15d，距下游邊界25d，當(dāng)重顆粒下降一個(gè)網(wǎng)格時(shí)，計(jì)算網(wǎng)格跟著重顆粒向下移動(dòng)一個(gè)網(wǎng)格。以上參數(shù)均在格子尺度下進(jìn)行設(shè)置。

2． 2 數(shù)值模型

采用單松弛格子Boltzmann方程求解流場(chǎng)，該模型是目前應(yīng)用最為廣泛的一種LB模型，其演化方程為

式中fi（x，t）為流體分子在i方向上的速度分布函數(shù)，fi（eq）（x，t）為對(duì)應(yīng)的平衡態(tài)分布函數(shù)，ci為i方向上的格子速度矢量，τ為無量綱的松弛時(shí)間因子。流體的宏觀速度u和密度ρf可以通過式（3）求解，

格子模型采用D2Q9格式，離散的速度矢量為

式中 wi為權(quán)系數(shù)，對(duì)于D2Q9為w0＝4／9，w1～4＝1／9，w5～8＝1／36。

對(duì)顆粒邊界的無滑移條件采用Lallemand等［13］提出的反彈格式來進(jìn)行處理，該方法可以看作是對(duì)標(biāo)準(zhǔn)反彈格式的一種改進(jìn)。另外，隨著顆粒的運(yùn)動(dòng)，一些固體節(jié)點(diǎn)從t到t＋Δt時(shí)刻會(huì)變?yōu)榱黧w節(jié)點(diǎn)，而一些流體節(jié)點(diǎn)則會(huì)變成固體節(jié)點(diǎn)，這對(duì)顆粒的運(yùn)動(dòng)會(huì)產(chǎn)生動(dòng)量變化，從而造成顆粒的受力變化，為了描述這部分力，本文采用Aidun等［14］提出的方法進(jìn)行處理。

3 驗(yàn) 證

采用通道中雙顆粒沉降的DKT結(jié)果驗(yàn)證本文計(jì)算程序的有效性。參數(shù)與文獻(xiàn)［15］相同，即顆粒直徑d＝0．2cm，顆粒密度為ρp＝1．01g／cm3，通道寬度L＝10d，高度為H＝40d，流體密度為ρf＝1g／cm3，粘性系數(shù)為ν＝10－6m2／s，初始時(shí)刻雙顆粒分別位于（5d，36d）和（5d，34d），分別記為Particle 1和Particle 2。計(jì)算采用的網(wǎng)格間距和時(shí)間步長分別為Δx＝10－4m和Δt＝5×10－4s，顆粒之間的碰撞采用了文獻(xiàn)［8，15］的彈性力模型。結(jié)果如圖1所示，雙顆粒在沉降過程中發(fā)生了典型的DKT現(xiàn)象，本文的計(jì)算結(jié)果與文獻(xiàn)［8，15］符合較好。

圖1 雙顆粒沉降的DKT結(jié)果對(duì)比Fig．1 Comparison of the DKT motion for sedimentation of two particles

4 計(jì)算結(jié)果及討論

針對(duì)5≤Re≤12范圍內(nèi)不同密度差異k的情況進(jìn)行數(shù)值研究。圖2給出了顆粒在Re＝5時(shí)的結(jié)果，圖中y1和y2分別為輕顆粒與重顆粒的橫向位置，且＝y(tǒng)1／d＝y(tǒng)2／d，t＊＝tU0／d。可以看出，當(dāng)k較小時(shí)，如k＝0或0．01，兩個(gè)顆粒在沉降中會(huì)發(fā)生周期性振蕩現(xiàn)象，且k越小，振幅越大。值得注意的是，k＝0意味著兩個(gè)顆粒完全相同，然而即使在初始位置對(duì)稱的情況下，這兩個(gè)顆粒也不會(huì)穩(wěn)定地沉降，此結(jié)果與文獻(xiàn)［2］符合?？梢钥闯?，當(dāng)k＝0．02時(shí)，兩個(gè)顆粒分別經(jīng)過初始階段的振蕩后逐漸穩(wěn)定并最終在通道中心沉降，即y＊1＝y(tǒng)＊2＝0。此外，計(jì)算結(jié)果表明Re＝6時(shí)的結(jié)果與圖2相似。

為了更好地對(duì)圖2所示的周期性運(yùn)動(dòng)進(jìn)行分析，本文通過兩個(gè)顆粒的水平位置坐標(biāo)y1和y2構(gòu)造出相關(guān)的相位圖，圖3分別給出了Re＝5和Re＝6時(shí)的結(jié)果?？梢钥闯?，每一個(gè)k對(duì)應(yīng)的相位圖均呈現(xiàn)出一條閉合曲線，意味著兩個(gè)顆粒做周期性運(yùn)動(dòng)。與圖5一致，隨著k的增大，閉合曲線尺寸逐漸減小。特別地，當(dāng)k＝0時(shí)，顆粒振蕩的軌跡相位圖呈現(xiàn)兩條曲線，如圖3（a）所示，這是典型的周期倍增現(xiàn)象。當(dāng)k值逐漸增大時(shí)，這種現(xiàn)象消失。此外，輕顆粒的振幅y1明顯大于重顆粒的振幅y2，而且輕顆粒幾乎均在通道左半邊沉降（y1＜0），而重顆粒則在通道右半邊沉降（y2＞0）。

圖2 Re＝5時(shí)顆粒水平位置的變化Fig．2 Horizontal displacements for the particles at Re＝5

從圖2可以看出，當(dāng)k較大時(shí)，顆粒會(huì)在通道中心穩(wěn)定地沉降，圖4給出了Re＝5時(shí)，0．02≤k≤0．08對(duì)應(yīng)的顆粒沉降流線圖?？梢钥闯觯?dāng)0．02≤k≤0．07時(shí)，兩個(gè)顆粒作為一個(gè)整體進(jìn)行沉降，如圖4（a～c）所示，這是由于輕顆粒處于重顆粒的尾流之中，使得阻力減小最終追上重顆粒。當(dāng)k足夠大時(shí)，如k＝0．08，重顆粒會(huì)離開輕顆粒，兩個(gè)顆粒最終獨(dú)立沉降，如圖4（d）所示。此外必須說明的是，圖中所示的箭頭表示顆粒的轉(zhuǎn)角，初始轉(zhuǎn)角為垂直方向。

圖3 顆粒運(yùn)動(dòng)軌跡的相位Fig．3 Trajectory of the particles in phase space

圖4 當(dāng)Re＝5時(shí)不同k值對(duì)應(yīng)的流線Fig．4 Instantaneous streamlines for different kat Re＝5

當(dāng)雷諾數(shù)增大時(shí)，顆粒的沉降形態(tài)會(huì)發(fā)生變化。圖5給出了Re＝7時(shí)顆粒的軌跡相位圖?？梢钥闯?，結(jié)果分成了兩個(gè)部分，k≤0．04時(shí)的結(jié)果與圖3類似；然而k繼續(xù)增大則出現(xiàn)了另一種沉降形態(tài)，即當(dāng)0．05≤k≤0．07時(shí)，輕顆粒在通道右半邊沉降（y1＞0），而重顆粒則正好相反（y2＜0），顯然，這兩個(gè)顆粒交換了初始位置，如圖5（b）所示。

當(dāng)k≥0．08時(shí)，顆粒仍然會(huì)穩(wěn)定地沉降。圖6給出了Re＝7時(shí)，k＝0．08，0．1和0．11對(duì)應(yīng)的流線圖。與圖4不同，當(dāng)k＝0．08和0．1時(shí)，兩個(gè)顆粒交換了初始位置，呈穩(wěn)定的錯(cuò)列結(jié)構(gòu)進(jìn)行沉降，如圖6（a，b）所示。此外，此時(shí)重顆粒轉(zhuǎn)動(dòng)的方向與輕顆粒一致，均是逆時(shí)針方向。造成這一現(xiàn)象的原因在于重顆?？拷诿娑艿捷^大的壁面排斥力，為了平衡這個(gè)力，導(dǎo)致了重顆粒的反常旋轉(zhuǎn)，并由此產(chǎn)生了一個(gè)馬格努斯力。

當(dāng)雷諾數(shù)繼續(xù)增大時(shí)，顆粒的沉降特性又會(huì)發(fā)生新的變化。如圖7所示，當(dāng)Re＝9時(shí)，顆粒沉降的軌跡相位圖也可以分成兩個(gè)部分，但與圖5所示的結(jié)果不同。當(dāng)k較大時(shí)，顆粒的振動(dòng)幅度隨k的增大而增大，而且明顯比k較小時(shí)的振動(dòng)幅度大（注意圖中坐標(biāo)軸的范圍），當(dāng)k＝0．09時(shí)觀察到了周期倍增的現(xiàn)象，如圖5（b）所示。此外，還可以看出，此時(shí)輕顆粒幾乎在通道右半邊沉降（y1＞0），而重顆粒則完全在右半邊沉降（y2＞0）。

圖5 Re＝7時(shí)顆粒運(yùn)動(dòng)軌跡的相位Fig．5 Trajectory of the particles in phase space at Re＝7

圖6 Re＝7時(shí)不同k值對(duì)應(yīng)的流線Fig．6 Instantaneous streamlines for different kat Re＝7

圖7 Re＝9時(shí)顆粒運(yùn)動(dòng)軌跡的相位Fig．7 Trajectory of the particles in phase space at Re＝9

為了更加直觀地觀察圖7對(duì)應(yīng)的兩類周期性沉降運(yùn)動(dòng)，圖8分別給出了0≤k≤0．05和0．06≤k≤0．08時(shí)的顆粒沉降軌跡。每個(gè)圖截取了兩個(gè)周期的顆粒位置與轉(zhuǎn)角，相鄰的時(shí)間間隔為500時(shí)間步?？梢钥闯觯?dāng)k≤0．05時(shí)，兩個(gè)顆粒在沉降中左右輕微振動(dòng)，相對(duì)位置沒有明顯變化，如圖8（a～c）所示。而當(dāng)k≥0．06時(shí)，情況發(fā)生本質(zhì)變化，重顆粒離開壁面并與輕顆粒交織在一起進(jìn)行沉降，如圖8（d～f）所示。當(dāng)重顆粒從壁面向通道中心運(yùn)動(dòng)時(shí)，輕顆粒處于其下游因此受到尾流的影響；當(dāng)輕顆粒加速逐漸靠近重顆粒時(shí)，二者之間的排斥作用逐漸增強(qiáng)，當(dāng)這個(gè)排斥力足夠大時(shí)，兩個(gè)顆粒又逐漸遠(yuǎn)離；而當(dāng)重顆?？拷诿鏁r(shí)，受到其排斥作用，又向中心運(yùn)動(dòng)，如此往復(fù)循環(huán)，形成了圖中所示的周期性運(yùn)動(dòng)。這種運(yùn)動(dòng)類似于典型的DKT過程，不同的是，由于這兩個(gè)顆粒密度差異較大，因此顆粒之間沒有發(fā)生翻轉(zhuǎn)。很明顯，當(dāng)k從0．05增大到0．06時(shí)，流體的慣性影響發(fā)生了本質(zhì)的變化，使得顆粒之間的相互作用發(fā)生了變化，造成了沉降特性的突變，如圖8（c，d）所示。當(dāng)進(jìn)一步增大雷諾數(shù)即當(dāng)Re＞9時(shí)，圖7（a）所示的沉降特性消失，只剩下類似于圖7（b）的結(jié)果，而從圖3和圖5可以看出，當(dāng)Re＜9時(shí)，本文只能觀察到與圖7（a）相似的結(jié)果。由此，可以將Re＝9看作是該沉降系統(tǒng)的一個(gè)臨界雷諾數(shù)。

圖9分別給出了Re＝10和Re＝12時(shí)對(duì)應(yīng)的顆粒沉降相位圖。對(duì)于Re＝10，僅顯示了k≥0．05的結(jié)果，這是因?yàn)樵谟?jì)算中發(fā)現(xiàn)當(dāng)k≤0．04時(shí)，兩個(gè)顆粒最終會(huì)穩(wěn)定地沉降，沉降的形態(tài)與圖6（a，b）相似，只是兩個(gè)顆粒的位置正好相反。從圖9可以看出，顆粒振動(dòng)的幅度隨k的增大而增大，尤其是重顆粒的幅度y2變化明顯，這與圖7（b）的結(jié)果相同，但此時(shí)相位圖出現(xiàn)了類似雙鈕線的形狀，這是因?yàn)榱黧w的慣性作用越來越顯著，使得顆粒的運(yùn)動(dòng)變得越來越復(fù)雜。

從上述研究可以看出，在沉降中一定條件下重顆粒與輕顆粒會(huì)一起沉降。圖10給出了不同雷諾數(shù)對(duì)應(yīng)的臨界k值即kc，只有當(dāng)兩個(gè)顆粒的密度差異大于該臨界值時(shí)，重顆粒才會(huì)擺脫輕顆粒并最終分離?？梢钥闯觯S著Re的增大，kc的值先增大后減小，當(dāng)Re在9附近時(shí)kc的值最大，這與上述研究中關(guān)于Re＝9為臨界雷諾數(shù)的結(jié)論一致。

圖8 Re＝9時(shí)顆粒沉降的軌跡Fig．8 Settling trajectory of particles at Re＝9

圖9 顆粒運(yùn)動(dòng)軌跡的相位Fig．9 Trajectory of the particles in phase space for different kat Re＝10and Re＝12

圖10 臨界值kc隨雷諾數(shù)的變化Fig．10 Dependence of the critical value kcon Re

5 結(jié) 論

本文采用基于動(dòng)量交換的格子Boltzmann方法對(duì)雙顆粒沉降系統(tǒng)進(jìn)行了直接數(shù)值模擬，研究了雷諾數(shù)Re＝5～12范圍內(nèi)顆粒之間密度差異k對(duì)顆粒沉降形態(tài)的影響，有以下結(jié)論。（1）當(dāng)雷諾數(shù)較小時(shí)，即當(dāng)5≤Re≤8時(shí)，顆粒之間的密度差k越小，顆粒周期性振動(dòng)的幅度越大，當(dāng)Re＝7時(shí)，兩個(gè)顆粒會(huì)交換初始位置繼續(xù)做周期性沉降運(yùn)動(dòng)。（2）當(dāng)雷諾數(shù)較大時(shí)，即當(dāng)10≤Re≤12時(shí)，由于流體的慣性作用，顆粒的振動(dòng)幅度顯著增大，而且k越大振幅越大；此外，計(jì)算發(fā)現(xiàn)存在一個(gè)臨界雷諾數(shù)，此時(shí)顆粒的沉降形態(tài)同時(shí)具有以上兩種特征。（3）當(dāng)雷諾數(shù)較大時(shí)，如果k較小，則兩個(gè)顆粒會(huì)呈現(xiàn)穩(wěn)定的錯(cuò)列結(jié)構(gòu)，此時(shí)重顆粒與輕顆粒的轉(zhuǎn)動(dòng)方向一致。此外，只有當(dāng)k大于某臨界值時(shí)，重顆粒才離開輕顆粒，此臨界值隨著Re先增大后減小。