劉 依， 聶德明

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

1 引 言

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

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

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

2 模型

2．1 物理模型

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

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

2． 2 數(shù)值模型

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

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

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

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

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

3 驗 證

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

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

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

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

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

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

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

圖3 顆粒運動軌跡的相位Fig．3 Trajectory of the particles in phase space

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

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

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

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

圖5 Re＝7時顆粒運動軌跡的相位Fig．5 Trajectory of the particles in phase space at Re＝7

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

圖7 Re＝9時顆粒運動軌跡的相位Fig．7 Trajectory of the particles in phase space at Re＝9

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

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

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

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

圖9 顆粒運動軌跡的相位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é) 論

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