丁 周，周 錕，孫 科，賀 鑄

(武漢科技大學(xué)省部共建耐火材料與冶金國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢，430081)

顆粒在流體中的自由沉降屬于固液兩相流[1-2]，普遍存在于自然界中，如泥沙的沉積[3-4]、大氣中水蒸氣以及霧霾顆粒的運(yùn)動(dòng)等[5-6]，同時(shí)，顆粒沉降還具有廣泛的工業(yè)生產(chǎn)背景，其相關(guān)研究在管道輸運(yùn)[7]、冶金化工流化床設(shè)備[8]等一些工程實(shí)際問題中也具有重要意義。顆粒沉降過程復(fù)雜，受到顆粒的粒徑、形狀、密度以及流體黏度等眾多因素影響。

針對一般的顆粒流動(dòng)問題，常用的數(shù)值模擬方法有歐拉-歐拉法(Euler-Euler method)[9]、歐拉-拉格朗日法(Euler-Lagrange method)[10]和顆粒解析法(particle-resolved method)[11]等三大類，這三類方法依次適用于顆粒非常小(如微納尺寸的氣溶膠顆粒)、顆粒比較小(如粉塵)以及顆粒相對比較大的情況。Peskin在研究心臟血管流動(dòng)時(shí)提出的浸入邊界法(immersed boundary method，IBM)[12-13]就是一種非常典型的顆粒解析法，該法的主要特點(diǎn)是利用統(tǒng)一、結(jié)構(gòu)化的歐拉網(wǎng)格來離散整個(gè)流固區(qū)域以及采用拉格朗日網(wǎng)格跟蹤流固界面，并借助delta函數(shù)實(shí)現(xiàn)二者的交互。傳統(tǒng)浸入邊界法基于彈性邊界，不適合對剛性邊界的顆粒流動(dòng)現(xiàn)象進(jìn)行數(shù)值模擬，而Uhlmann[14]則通過引入直接力所提出的直接力-浸入邊界法(direct forcing-immersed boundary method, DF-IBM)很好地解決了這一問題，借助DF-IBM對二維條件下的雙圓形顆粒沉降以及三維條件下單個(gè)和多個(gè)球形顆粒沉降進(jìn)行研究時(shí)，模擬結(jié)果與相關(guān)的實(shí)驗(yàn)結(jié)果[15]基本吻合，Di等[16]采用DF-IBM對經(jīng)典的二維雙圓形顆粒的牽引、碰撞及翻滾(DKT)[17]進(jìn)行模擬，也獲得了較理想的研究結(jié)果，此外，還有宮兆新等[18]對浸入邊界法進(jìn)行了較全面的研究。一直以來，在DF-IBM中歐拉網(wǎng)格和拉格朗日網(wǎng)格的尺寸被認(rèn)為應(yīng)該相互匹配，并且對于每個(gè)拉格朗日節(jié)點(diǎn)都需要定義一個(gè)與歐拉網(wǎng)格匹配的權(quán)重積分，然而本課題組最新的研究[19]表明，拉格朗日體積權(quán)重?zé)o需特殊選擇，可以在很大范圍內(nèi)自由取值。

當(dāng)前，有關(guān)二維條件下雙圓形顆粒和三維條件下雙球顆粒的沉降實(shí)驗(yàn)與模擬方法已經(jīng)比較完備，但針對多顆粒的沉降問題仍停留在宏觀分析階段，考慮到實(shí)際顆粒輸運(yùn)過程中往往是多個(gè)顆粒一起沉降或者提升，故本文分別以單個(gè)球形、線形排列的雙球和多球顆粒為研究對象，借助直接力浸入邊界法對相應(yīng)的沉降過程進(jìn)行數(shù)值模擬研究，并對剛體假設(shè)(rigid body assumption,RBA)[20]的結(jié)果準(zhǔn)確性進(jìn)行分析，重點(diǎn)探討了線形排列多球顆粒的沉降行為。

1 模型與方法

1.1 物理模型

(a)模擬區(qū)域 (b) 球形顆粒表面離散點(diǎn)圖1 物理模型Fig.1 Physical model

1.2 直接力浸入邊界法

采用直接力浸入邊界法對顆粒自由沉降進(jìn)行數(shù)值模擬，其數(shù)理模型中不可壓縮流體的控制方程為

(1)

(2)

式(1)～式(2)中:u、p分別為流體的速度和壓力；為哈密頓算子；fIBM為數(shù)值模擬中歐拉網(wǎng)格點(diǎn)處的作用力。顆粒在沉降時(shí)的運(yùn)動(dòng)方程為

(3)

式中：mP為顆粒質(zhì)量；t為時(shí)間；τ為應(yīng)力張量，τ=-Ip+vf(u+uT)，其中I是二階單位張量,T表示矢量矩陣的轉(zhuǎn)置；n為顆粒表面法線方向；GP為顆粒所受重力，GP=VP(ρf-ρP)，其中VP為顆粒體積，ρP為顆粒密度。式(3)中的相互作用力項(xiàng)

(4)

式中，χ表示歐拉節(jié)點(diǎn)。采用RBA近似處理式(4)等號(hào)右邊第2項(xiàng)，有

(5)

(6)

式中，fIBM可表示為

fIBM=BWFIBM

(7)

式中:B為插值函數(shù)，W為控制收斂的權(quán)重因子，F(xiàn)IBM為顆粒表面拉格朗日離散點(diǎn)上的作用力，也就是說fIBM可由顆粒表面拉格朗日離散點(diǎn)上的作用力FIBM經(jīng)插值處理而獲得，而FIBM又可由顆粒與流體在界面上的速度差與時(shí)間步長的比率來確定，即

BiI=ψ(χi-ΧI)

(8)

(9)

(10)

(11)

(12)

需要指出的是，本研究模擬顆粒尺寸遠(yuǎn)大于流體分子的平均自由程，因此顆粒表面所受流體分子的隨機(jī)碰撞處于平衡狀態(tài)，也就是說，流體的分子熱運(yùn)動(dòng)對顆粒的影響是各向均勻的，其宏觀作用表現(xiàn)為流體的黏性。對于液體來說，其分子平均自由程為納米量級，而對于常溫下的空氣，其分子平均自由程為亞微米量級，所以，本文模擬適用于顆粒尺寸為微米量級以上的情況，文中顆粒尺寸單位均為其尺寸量級所對應(yīng)的國際單位制單位及其導(dǎo)出單位。

2 模擬方法驗(yàn)證

為進(jìn)一步驗(yàn)證數(shù)值模擬對網(wǎng)格的依賴性，在3組不同網(wǎng)格密度條件下對類似的顆粒沉降問題進(jìn)行了模擬，模擬結(jié)果見圖3。其中NP表示顆粒拉格朗日離散點(diǎn)數(shù)目，不同的拉格朗日離散點(diǎn)數(shù)對應(yīng)不同密度尺寸的歐拉網(wǎng)格，相應(yīng)的誤差分析如表2所示。結(jié)合圖3及表2結(jié)果可知，當(dāng)顆粒拉格朗日離散點(diǎn)為258，相應(yīng)槽道歐拉網(wǎng)格密度為68×544×69即Dp/h=9.063時(shí)，模擬計(jì)算結(jié)果較其它網(wǎng)格密度條件下的模擬值更為準(zhǔn)確，采取此網(wǎng)格進(jìn)行模擬計(jì)算，模擬過程中Courant-Friedrichs-Lewy(CFL)數(shù)均小于1。

表1 算例的對比結(jié)果

(a)算例1

(b)算例2

(c)算例3

圖3 不同網(wǎng)格密度條件下的顆粒沉降速度Fig.3 Particle sedimentation velocities under various grid densities

表2 誤差分析

3 結(jié)果及討論

3.1 剛體假設(shè)(RBA)的影響

基于文獻(xiàn)[14]中的算例2，通過計(jì)算，分析了不同RBA條件下，式(5)中積分項(xiàng)(等號(hào)左端)及乘積項(xiàng)(等號(hào)右端)的數(shù)值隨時(shí)間變化的規(guī)律，結(jié)果如圖4所示。圖4中，標(biāo)記為RBA乘積值和非RBA乘積值的曲線分別反映了顆粒是否為剛體時(shí)顆粒體積與顆粒質(zhì)心沉降速度的乘積隨時(shí)間變化的趨勢，標(biāo)記為RBA積分值和非RBA積分值的曲線分別反映了顆粒是否為剛體時(shí)借助Monte Carlo隨機(jī)積分法計(jì)算所得式(5)中積分項(xiàng)的值隨時(shí)間變化的趨勢，作為參照，標(biāo)記為RBA實(shí)驗(yàn)乘積值的曲線則反映了顆粒為剛體時(shí)，由文獻(xiàn)[15]實(shí)驗(yàn)值所得出的顆粒質(zhì)心沉降速度與顆粒體積的乘積隨時(shí)間變化的趨勢。上述積分項(xiàng)、乘積項(xiàng)值隨時(shí)間的變化趨勢即為顆粒沉降速度隨時(shí)間的變化趨勢。由圖4可見，將顆?？醋鲃傮w時(shí)，相應(yīng)的計(jì)算結(jié)果更接近于文獻(xiàn)[15]所報(bào)道的實(shí)驗(yàn)值；當(dāng)顆粒不是剛體時(shí)，計(jì)算值和實(shí)驗(yàn)值的誤差較大。此外，將顆粒作為剛體處理時(shí)，相應(yīng)的乘積值曲線相較積分值曲線更為平滑、波動(dòng)更小。以上結(jié)果表明，在球狀顆粒沉降過程中，當(dāng)顆粒為剛體時(shí)，相應(yīng)計(jì)算結(jié)果更加準(zhǔn)確。

圖4 不同RBA條件下方程(5)的計(jì)算值Fig.4 Calculated values of equation (5) under different RBA conditions

3.2 雙球顆粒沉降

在雙球顆粒沉降模擬中，設(shè)定2個(gè)球體顆粒間球心的初始間距分別為1DP、1.5DP和2DP。在顆粒沉降過程中，不同初始間距條件下，2顆粒相對間距Dr隨時(shí)間t變化的模擬結(jié)果如圖5所示。由圖5可見，在球狀顆粒沉降過程初始階段，2球間發(fā)生相互牽引致使二者間距逐漸縮短，這與經(jīng)典雙球沉降模型[17]中的牽引過程一致，當(dāng)2球相對距離減小到一定數(shù)值時(shí)則會(huì)發(fā)生2球相互觸碰，初始間距不同的2球均能發(fā)生牽引和觸碰，且發(fā)生觸碰的時(shí)刻隨2球初始間距的增大而推遲。隨著顆粒沉降的繼續(xù)進(jìn)行，2球相互觸碰后又發(fā)生翻滾運(yùn)動(dòng)導(dǎo)致二者分離，2球相對間距又開始逐漸增大(圖5(a))。

(a)初始間距1DP

(b)初始間距1.5DP

(c)初始間距2DP

圖6所示為雙球顆粒沉降過程中的速度云圖及壓力云圖。需要特別指明的是，在DF-IBM方法中，由于采用統(tǒng)一的歐拉網(wǎng)格來處理整個(gè)區(qū)域(包含顆粒所占據(jù)的空間)，因此也有流體填充在顆粒內(nèi)部，而顆粒內(nèi)部的速度場，并不完全滿足顆粒的剛體運(yùn)動(dòng)(見圖6(a))，此問題的根源還有待進(jìn)一步深入研究。從雙球的壓力云圖(圖6(b))中可以清楚看出，單個(gè)球體的下方為正壓區(qū)，上方為負(fù)壓區(qū)。當(dāng)2號(hào)球處于1號(hào)球的負(fù)壓區(qū)時(shí)，2號(hào)球所受壓差阻力減小并導(dǎo)致其加速沉降，進(jìn)一步的分析表明，1號(hào)球的沉降速度較其單獨(dú)沉降時(shí)要小，也就是說，2號(hào)球在1號(hào)球的牽引下加速運(yùn)動(dòng)，最終追上后者發(fā)生觸碰。

(a)速度 (b)壓力圖6 雙球顆粒沉降云圖Fig.6 Sedimentation nephograms of double spherical particles

3.3 線形排列多球顆粒沉降

線形排列多球模型之于雙球模型，可看作在雙球上方等間距放入一定數(shù)量的小球，原有的雙球模型在加入更多球形顆粒后成為具有多球系統(tǒng)的新模型。多球顆粒系統(tǒng)沉降過程中的速度云圖及壓力云圖如圖7所示。從圖7中可以看出，當(dāng)顆粒已經(jīng)沉降一段距離后，在球形顆粒數(shù)分別為3、4、5、6的算例中，1號(hào)球和2球的間距均小于其它相鄰兩球之間的間距，這表明1號(hào)球和2號(hào)球率先發(fā)生牽引，并且在不同的球形顆粒數(shù)條件下，1號(hào)球下方的正壓梯度均大于其它球的相應(yīng)值。相比雙球模型，多球系統(tǒng)中的2號(hào)球上方還存在由3號(hào)球造成的正壓區(qū)，因此，在3號(hào)球下方正壓區(qū)及1號(hào)球上方負(fù)壓區(qū)的共同作用下，2號(hào)球會(huì)更快地追上1號(hào)球。

(a)速度

(b)壓力

因?yàn)楫?dāng)線形排列的球狀顆粒數(shù)超過3以后，伴隨著1號(hào)、2號(hào)球發(fā)生DKT(牽引，觸碰、翻滾)過程，1～3號(hào)球的運(yùn)動(dòng)狀態(tài)基本一致，故而對1～3號(hào)球中相鄰球間的相對間距變化進(jìn)行分析，結(jié)果如圖8所示，其中P2-P1表示1號(hào)球和2號(hào)球的相對間距，P3-P2表示2號(hào)球和3號(hào)球的相對間距。由圖8可見，在顆粒沉降過程中，1號(hào)球和2號(hào)球先發(fā)生牽引和觸碰，然后再分離，這與二者在雙球模型中的表現(xiàn)一致，區(qū)別在于整個(gè)過程因3號(hào)球的加入而導(dǎo)致進(jìn)度加快。同時(shí)，3號(hào)球和2號(hào)球僅僅發(fā)生牽引而不會(huì)產(chǎn)生觸碰。此外，本研究模擬結(jié)果還表明，隨著線形排列顆粒數(shù)目的不斷增加，使用DF-IBM對相應(yīng)沉降過程進(jìn)行模擬計(jì)算的效率始終保持在較高水平。IBM模擬主要采用歐拉網(wǎng)格及拉格朗日網(wǎng)格，其中歐拉網(wǎng)格的位置和數(shù)量是固定不變的，而拉格朗日網(wǎng)格的變化則依賴于顆粒數(shù)目。在模擬計(jì)算中，計(jì)算量主要取決于歐拉網(wǎng)格的數(shù)量，這是因?yàn)樵谡麄€(gè)使用了壓力投影法的模擬過程中，對壓力場的泊松方程進(jìn)行求解耗時(shí)最多，而泊松方程又是在歐拉網(wǎng)格上離散求解的。一般來說，顆粒越多，計(jì)算量越大，但是IBM法中歐拉網(wǎng)格數(shù)量固定，所以計(jì)算量的增加幅度非常有限，這也是IBM方法目前應(yīng)用比較廣泛的原因之一。

圖8 多球顆粒模型的Dr-t曲線

3.4 線形排列雙球與多球顆粒沉降對比

借助雙球和三球顆粒模型對比分析了沉降過程中1號(hào)球和2號(hào)球的相對間距與碰撞時(shí)間的關(guān)系，結(jié)果如圖9所示。由圖9可見，當(dāng)相鄰顆粒之間的初始間距從1DP增大到4DP時(shí)，即隨著Dr的增加，雙球和三球模型中的1號(hào)球和2號(hào)球均能發(fā)生觸碰，而且在三球模型中，因3號(hào)球的引入導(dǎo)致1號(hào)、2號(hào)球觸碰時(shí)間提前。

當(dāng)相鄰兩球顆粒初始間距為1DP時(shí)，在雙球和三球模型中發(fā)生牽引的1號(hào)球和2號(hào)球的沉降速度隨二者相對間距變化的曲線如圖10所示，其中2P-P1、3P-P1分別表示雙球和三球模型中的1號(hào)球，即牽引球；2P-P2和3P-P2分別表示雙球和三球模型中的2號(hào)球，即受牽引球。從圖10中可以看出，在沉降開始后的一段時(shí)間內(nèi)，兩種模型中1號(hào)球和2號(hào)球都在加速沉降且二者相對間距基本不變，緊接著雙球模型中的兩球間距稍有增大然后再減小，而與此同時(shí)，三球模型中發(fā)生牽引的1號(hào)球和2號(hào)球間的相對間距卻一直在縮小，表明三球模型中的牽引現(xiàn)象較雙球模型先發(fā)生，應(yīng)歸因于三球模型中3號(hào)球的引入加速了這一過程。

圖9 雙球和三球顆粒模型的碰撞時(shí)間Fig.9 Collision time of double and three spherical particles models

圖10 雙球和三球顆粒模型的沉降速度Fig.10 Sedimentation velocities of double and three spherical particles models

在球形顆粒沉降過程中，針對1號(hào)球和2號(hào)球在三球模型中較雙球模型先發(fā)生牽引的問題，從球形顆粒所受阻力角度進(jìn)行分析，兩種模型中1號(hào)球及2號(hào)球所受阻力變化情況如圖11所示。由圖5與圖8結(jié)果可知，沉降時(shí)間t超過0.5后，1號(hào)球和2號(hào)球之間的相對間距出現(xiàn)明顯變化，而從圖11(a)中又可以看出，受牽引的2號(hào)球在三球模型中所受阻力較其在雙球模型中所受阻力要小，從而能在更短的時(shí)間內(nèi)追上1號(hào)牽引球以實(shí)現(xiàn)牽引。在顆粒沉降初始階段，1號(hào)牽引球在三球模型中所受阻力與其在雙球模型中所受阻力幾乎一樣，當(dāng)t超過0.7時(shí)，1號(hào)牽引球在三球模型及雙球模型中所受阻力出現(xiàn)明顯差別(圖11(b))，這是因?yàn)橄嗤瑫r(shí)刻三球模型中1號(hào)球和2號(hào)球相對間距小于雙球模型中的相應(yīng)值，導(dǎo)致1號(hào)球所受斥力增大而阻力減小。

(a)受牽引球(2號(hào)球)

(b)牽引球(1號(hào)球)

4 結(jié)語

當(dāng)球形顆粒初始間距介于一定范圍內(nèi)且視球形顆粒為剛體時(shí)，豎直線形排列的雙球或多球顆粒(顆粒按位置從低到高依次排序)中的1、2號(hào)球均可進(jìn)行DKT(牽引、觸碰、翻滾)運(yùn)動(dòng)。在雙球顆?；A(chǔ)上引入數(shù)量不等的球形顆粒形成線形排列多球顆粒，能改變2號(hào)球的受力從而加速1號(hào)球和2號(hào)球DKT運(yùn)動(dòng)的發(fā)生，這一現(xiàn)象對引入顆粒的數(shù)量不敏感。本研究方法可直接擴(kuò)展應(yīng)用到任意非規(guī)則形狀顆粒的沉降研究，能為顆粒沉降問題分析提供一定的參考。