葉 挺,凡鳳仙

（上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093）

管道內(nèi)顆粒負(fù)載流動(dòng)是能源、化工及生命科學(xué)等領(lǐng)域的研究熱點(diǎn)之一。顆粒橫向遷移是管道內(nèi)顆粒負(fù)載流動(dòng)中的一種有趣現(xiàn)象，它是指隨機(jī)散布在泊肅葉流中的中性浮力顆粒（密度與流體密度相同的顆粒），經(jīng)過一定時(shí)間后，遷移到距管中心一定距離橫向位置上流動(dòng)的現(xiàn)象[1]。顆粒橫向遷移現(xiàn)象有望在顆粒的分離、篩選、混合、反應(yīng)中得到應(yīng)用[2-3]，富有研究價(jià)值和應(yīng)用潛力。顆粒橫向遷移的早期研究主要采用實(shí)驗(yàn)方法，且多關(guān)注雷諾數(shù)、顆粒尺寸等對(duì)顆粒遷移的最終平衡位置的影響。隨著計(jì)算機(jī)和數(shù)值模擬技術(shù)的發(fā)展，數(shù)值模擬逐漸成為研究顆粒負(fù)載流動(dòng)的重要手段，通過數(shù)值模擬可以方便地對(duì)參數(shù)進(jìn)行調(diào)整，且能夠獲取豐富的流場和顆粒動(dòng)態(tài)過程的細(xì)節(jié)信息，便于揭示相關(guān)現(xiàn)象的規(guī)律及內(nèi)在機(jī)理。近年來，多種數(shù)值模擬方法在顆粒橫向遷移研究中得到應(yīng)用，并取得了一些有價(jià)值的研究成果，但由于顆粒橫向遷移所涉及的液固兩相流場的復(fù)雜性，目前對(duì)顆粒橫向遷移行為的認(rèn)識(shí)仍存在爭議，對(duì)橫向遷移內(nèi)在機(jī)理研究還不夠深入。本文首先針對(duì)泊肅葉流下中性浮力顆粒橫向遷移的相關(guān)實(shí)驗(yàn)和數(shù)值模擬研究進(jìn)行歸納和評(píng)述，然后基于直接力/虛擬區(qū)域（direct-forcing fictitious domain，DF/FD）方法模擬再現(xiàn)單個(gè)顆粒橫向遷移的過程，在此基礎(chǔ)上指出今后進(jìn)行顆粒橫向遷移研究的思路和方向。

1 泊肅葉流下中性浮力顆粒橫向遷移的實(shí)驗(yàn)研究

1.1 圓管內(nèi)顆粒橫向遷移的實(shí)驗(yàn)研究

1961年，Segré和Silberberg[4]首次在實(shí)驗(yàn)中觀察到圓管內(nèi)低雷諾數(shù)（7

圖1 橫向遷移后顆粒濃度沿徑向的分布Fig. 1 Particle concentration as a function of radial position after lateral migration

隨后，一些學(xué)者迅速開展了不同條件下圓管內(nèi)顆粒橫向遷移平衡位置的探討。1962年，Oliver[5]通過實(shí)驗(yàn)發(fā)現(xiàn)顆粒的旋轉(zhuǎn)對(duì)于顆粒遷移的最終穩(wěn)定位置有顯著影響，無旋轉(zhuǎn)時(shí)顆粒最終穩(wěn)定位置更接近管中心軸線，并認(rèn)為向內(nèi)的壁面效應(yīng)與向外的馬格努斯效應(yīng)使得顆粒最終穩(wěn)定在特定的徑向位置。1965年，Jeffrey等[6]在圓管直徑D=32.5 mm，顆粒直徑d=1.5，2.0，2.9 mm的條件下，實(shí)驗(yàn)觀察到中性浮力顆粒最終遷移到距管中心軸線約為2R/3的環(huán)形區(qū)域。1966年，Karnis等[7]通過實(shí)驗(yàn)發(fā)現(xiàn)中性浮力剛性顆粒橫向遷移的平衡位置與顆粒形狀、粒徑和管徑比有關(guān)，中性浮力可變形顆粒總是遷移到管中心軸線位置。

然而，由于受實(shí)驗(yàn)技術(shù)的限制，難以得到更為詳細(xì)的顆粒動(dòng)力學(xué)信息和流場信息，相關(guān)實(shí)驗(yàn)研究一度停滯。直到進(jìn)入21世紀(jì)，微流控技術(shù)興起，泊肅葉流下中性浮力顆粒的橫向遷移重新引起了學(xué)者的重視。2004年，Eloot等[8]實(shí)驗(yàn)研究了顆粒直徑d=10 μm的近中性浮力球形顆粒在直徑D=220，530 μm的毛細(xì)管內(nèi)的運(yùn)動(dòng)，結(jié)果表明，在較高的雷諾數(shù)、較長和較細(xì)的毛細(xì)管、較稀的顆粒濃度的條件下，顆粒橫向遷移的平衡位置更接近于管壁。2004年，Matas等[1]實(shí)驗(yàn)研究了雷諾數(shù)Re在67～1700范圍以及達(dá)到2 000時(shí)，圓管內(nèi)稀相泊肅葉流下中性浮力球形顆粒（D/d=8～42）的橫向遷移，發(fā)現(xiàn)Re較小時(shí)，顆粒遷移到距管中心軸線一定距離的環(huán)形區(qū)域，隨著Re的增大，環(huán)形區(qū)域的位置向管壁移動(dòng)，當(dāng)Re繼續(xù)增大到一定值時(shí)，顆粒的遷移使得顆粒分布在2個(gè)環(huán)形區(qū)域，即除了靠近管壁的外環(huán)外，還存在一個(gè)靠近管中心的內(nèi)環(huán)。圖2給出了Matas等[1]實(shí)驗(yàn)得到的Re=1000時(shí)顆粒概率密度P沿徑向的分布情況，可見內(nèi)環(huán)出現(xiàn)在r/R≈0.5處，外環(huán)出現(xiàn)在r/R≈0.9處。2017年，Morita等[9]實(shí)驗(yàn)研究了雷諾數(shù)Re在100～1100范圍、顆粒粒徑與管徑比約為0.1時(shí)，泊肅葉流下中性浮力顆粒橫向遷移的平衡位置，研究發(fā)現(xiàn)，Re較大時(shí)，顆粒分布首先集中在2個(gè)環(huán)形（即內(nèi)環(huán)和外環(huán)）區(qū)域，隨著顆粒向下游運(yùn)動(dòng)，內(nèi)環(huán)附近顆粒減少，而外環(huán)附近顆粒增加，若管道足夠長，顆粒將最終遷移到外環(huán)位置。Morita等[9]得到了較大雷諾數(shù)下顆粒遷移特性的演變過程示意圖，如圖3所示，即隨機(jī)分布的顆粒隨流體進(jìn)入管內(nèi)（x=0），首先向內(nèi)遷移到達(dá)一定的位置（x=L2）；接著轉(zhuǎn)而向外遷移，引起顆粒在內(nèi)環(huán)和外環(huán)2個(gè)區(qū)域集中（x=L1）；顆粒繼續(xù)向外遷移，最終集中在一個(gè)環(huán)形區(qū)域（x=L0）。在圖3中，顆粒向內(nèi)遷移發(fā)生在流體從初始狀態(tài)充分發(fā)展為泊肅葉流的管段，這表明顆粒起初向內(nèi)遷移與流動(dòng)的發(fā)展過程有關(guān)，顆粒轉(zhuǎn)而向外遷移則是由于泊肅葉流中橫向速度梯度引起的升力和顆粒的旋轉(zhuǎn)引起的升力所致。

圖2 橫向遷移后顆粒概率密度沿徑向的分布Fig.2 Particle probability density as a function of radial position after lateral migration

圖3 較大雷諾數(shù)下顆粒遷移特性的演變過程示意圖Fig.3 Schematic diagram of the evolution of particle migration at relatively high Reynolds number

1.2 矩形與方形通道內(nèi)的實(shí)驗(yàn)研究

由于壁面效應(yīng)的影響，非中心對(duì)稱的矩形和方形截面通道內(nèi)顆粒橫向遷移的特性將與中心對(duì)稱的圓管不同，一些學(xué)者針對(duì)矩形與方形截面通道內(nèi)顆粒的橫向遷移開展了一系列研究。

在矩形截面通道泊肅葉流下顆粒橫向遷移的研究方面，2005年，Staben等[10]實(shí)驗(yàn)研究了雷諾數(shù)Re（以矩形窄邊長H作為特征尺寸）在0.005～0.05范圍時(shí)顆粒粒徑對(duì)不同管口形狀的矩形截面微通道內(nèi)顆粒橫向遷移的影響，結(jié)果表明，小粒徑顆粒（0.07≤d/H≤0.10）在直管口通道內(nèi)向通道中心遷移，在帶倒角管口通道內(nèi)分布更為均勻；中等尺度粒徑顆粒（0.46≤d/H≤0.54）在直管口和帶倒角管口通道內(nèi)均呈幾乎均勻的分布；少量大粒徑顆粒（0.79≤d/H≤0.95）的分布略微向通道的一個(gè)側(cè)壁傾斜。2008年，Bhagat等[11]在Re>1的條件下實(shí)驗(yàn)研究了矩形截面微通道內(nèi)顆粒的橫向遷移特性，結(jié)果表明，顆粒遷移過程主要取決于通道截面較短邊的幾何尺度，顆粒最終分布在較大的壁面附近，d/H=0.07情況下，Re≥ 10時(shí)才能發(fā)生橫向遷移，且隨著d/H的增加，顆粒橫向遷移所需要的臨界雷諾數(shù)降低。由于文獻(xiàn)[10-11]的實(shí)驗(yàn)中采用的雷諾數(shù)范圍有很大差異，使得實(shí)驗(yàn)結(jié)果不同，雷諾數(shù)對(duì)矩形微通道內(nèi)顆粒橫向遷移的影響仍需要進(jìn)一步系統(tǒng)的研究。

在方形截面通道泊肅葉流下顆粒橫向遷移的研究方面，2008年，Kim等[12]實(shí)驗(yàn)研究了通道雷諾數(shù)Re在0.06～58.65范圍內(nèi)中性浮力顆粒（通道水力直徑與粒徑之比λ≈14）的橫向遷移，結(jié)果表明，在很低的雷諾數(shù)下，顆粒的橫向遷移仍會(huì)發(fā)生，存在一個(gè)在20～30之間的臨界雷諾數(shù)Rec，當(dāng)ReRec時(shí)，顆粒遷移的平衡位置向管壁移動(dòng)。2009年，Cho等[13]實(shí)驗(yàn)研究了通道雷諾數(shù)Re在0.25～2.5范圍時(shí)中性浮力顆粒（λ≈5，8，13.3）的橫向遷移，結(jié)果表明，顆粒遷移受到雷諾數(shù)Re和顆粒粒徑的影響，顆粒粒徑過小（λ≈13.3）時(shí)，不能發(fā)生顆粒遷移。2014年，Miura等[14]實(shí)驗(yàn)研究了雷諾數(shù)Re在100～1200范圍時(shí)中性浮力顆粒（通道邊長與粒徑之比為9.23）的橫向遷移情況，研究發(fā)現(xiàn)：Re>250時(shí)顆粒最終穩(wěn)定在截面中心軸線和轉(zhuǎn)角的8個(gè)平衡位置上，且隨著雷諾數(shù)的增大，管截面轉(zhuǎn)角處的平衡位置將遠(yuǎn)離通道中心，而管截面中心軸線處的平衡位置則將靠近通道中心；Re<250時(shí)轉(zhuǎn)角處的平衡位置消失。典型雷諾數(shù)下顆粒遷移后的最終穩(wěn)定位置如圖4所示，圖4中x和y表示位置坐標(biāo)，a為通道橫截面邊長。可見，方形截面通道內(nèi)顆粒的橫向遷移特性與圓管存在顯著差異。

圖4 典型雷諾數(shù)下橫向遷移后的顆粒分布Fig.4 Particle distributions after lateral migration at typical Reynolds numbers

1.3 其他相關(guān)研究

為了理解顆粒橫向遷移的機(jī)理，1979年，Aokl等[15]實(shí)驗(yàn)研究了基于顆粒自由沉降速度的雷諾數(shù)ReF（ReF=2r'uF/ν，r'為顆粒半徑，uF為顆粒在靜止流體中的自由沉降速度，ν為流體運(yùn)動(dòng)黏度）在0.043～10范圍時(shí)非中性浮力顆粒在圓管內(nèi)的橫向遷移，結(jié)果表明，ReF>1時(shí)，若顆粒運(yùn)動(dòng)滯后于流體，則顆粒向管中心軸線遷移，若顆粒運(yùn)動(dòng)超前于流體，則顆粒向管壁遷移，這一現(xiàn)象受控于Magnus效應(yīng)；ReF<1時(shí)，靠近管壁的顆粒向內(nèi)遷移，靠近管中心的顆粒向外遷移，最終顆粒到達(dá)管壁和管中心軸線之間的一個(gè)平衡位置，這一現(xiàn)象是由壁面效應(yīng)和慣性作用共同引起的。為了探討顆粒橫向遷移在微流控領(lǐng)域的應(yīng)用，2011年，文獻(xiàn)[16]基于螺旋微通道內(nèi)慣性升力和Dean曳力共同作用下，不同粒徑顆粒發(fā)生橫向遷移的最終平衡位置不同的原理，設(shè)計(jì)出利用螺旋形通道富集和分離顆粒的微流控裝置。2012年，Martel等[17]實(shí)驗(yàn)研究了高寬比在0.125～1范圍的螺旋形微通道內(nèi)顆粒富集動(dòng)力學(xué)行為，給出了實(shí)現(xiàn)高質(zhì)量富集所需要的通道長度和流速范圍。這些研究可以為理解顆粒橫向遷移的機(jī)理以及開發(fā)基于橫向遷移原理的微流控裝置提供參考。

2 泊肅葉流下中性浮力顆粒橫向遷移的數(shù)值模擬研究

與實(shí)驗(yàn)研究相比，管道內(nèi)中性浮力顆粒橫向遷移的數(shù)值模擬研究開展得較晚，但由于數(shù)值模擬的迅速發(fā)展以及計(jì)算機(jī)能力的快速提升，數(shù)值模擬已成為研究顆粒橫向遷移行為及其內(nèi)在機(jī)理的重要手段。目前，研究顆粒橫向遷移的數(shù)值模擬方法主要有有限元方法、虛擬區(qū)域方法、格子波爾茲曼方法及界面追蹤方法等。

2.1 有限元方法

1994年，F(xiàn)eng等[18]對(duì)泊肅葉流下中性浮力圓形顆粒的運(yùn)動(dòng)進(jìn)行二維有限元模擬，獲得了與文獻(xiàn)[4]的實(shí)驗(yàn)類似的顆粒橫向遷移效應(yīng)，并將顆粒橫向遷移的驅(qū)動(dòng)力歸因于潤滑產(chǎn)生的壁面排斥力、剪切滑移產(chǎn)生的慣性升力、顆粒旋轉(zhuǎn)產(chǎn)生的升力和速度曲線的曲率產(chǎn)生的升力。

2009年，文獻(xiàn)[19]利用有限元方法求解穩(wěn)態(tài)納維?斯托克斯方程，通過表面積分的方法獲得顆粒受到的力和力矩，對(duì)狹窄矩形截面微通道內(nèi)顆粒橫向遷移進(jìn)行數(shù)值模擬，并進(jìn)行了相應(yīng)的實(shí)驗(yàn)，數(shù)值模擬得到的顆粒橫向遷移的平衡位置與實(shí)驗(yàn)結(jié)果吻合良好。在此基礎(chǔ)上，模擬了顆粒橫向遷移過程中的升力，分析了升力的標(biāo)度關(guān)系。

2.2 虛擬區(qū)域方法

2005年，Yang等[20]采用虛擬區(qū)域方法，通過固定網(wǎng)格分布式拉格朗日乘子描述顆粒的運(yùn)動(dòng)，模擬了管內(nèi)單個(gè)中性浮力剛性顆粒的橫向遷移，得到的顆粒平衡位置與實(shí)驗(yàn)吻合良好。2008年，Pan等[21]采用基于拉格朗日乘子的虛擬區(qū)域方法模擬三維泊肅葉流中的中性浮力橢圓形顆粒的運(yùn)動(dòng)，研究了顆粒的旋轉(zhuǎn)和取向特性，發(fā)現(xiàn)隨著雷諾數(shù)和顆粒形狀的變化，顆粒呈現(xiàn)截然不同的運(yùn)動(dòng)狀態(tài)，顆粒橫向遷移的平衡位置與文獻(xiàn)報(bào)道的圓盤形顆粒的實(shí)驗(yàn)結(jié)果[7]類似。上述方法又稱為DLM/FD（distributed-Lagrange-multiplier fictitious domain）方法，該方法需要顯式計(jì)算拉格朗日乘子進(jìn)而計(jì)算顆粒的運(yùn)動(dòng)參數(shù)，因而計(jì)算效率不高。

為了提高虛擬區(qū)域方法的計(jì)算效率，直接力/虛擬區(qū)域（direct-forcing fictitious domain，DF/FD）方法被提出，該方法引入體積力以計(jì)算顆粒的運(yùn)動(dòng)參數(shù)，從而避免了對(duì)拉格朗日乘子的依賴。2008年，Shao等[22]在流動(dòng)方向上施加周期性邊界條件，利用DF/FD方法模擬了圓管泊肅葉流下球形中性浮力顆粒的運(yùn)動(dòng)過程以及顆粒運(yùn)動(dòng)對(duì)流場的影響，發(fā)現(xiàn)在雷諾數(shù)高于臨界值時(shí)顆粒橫向遷移到一定的平衡位置，這與文獻(xiàn)[1]的實(shí)驗(yàn)結(jié)果在定性上一致。2009年，Sun等[23]采用DF/FD方法模擬了二維通道內(nèi)的圓形中性浮力顆粒在恒定壓力和波動(dòng)壓力驅(qū)動(dòng)的流動(dòng)中，橫向遷移過程以及顆粒附近的流場狀況，但缺少和實(shí)驗(yàn)結(jié)果的對(duì)比驗(yàn)證。

2.3 格子玻爾茲曼方法

2006年，Chun等[24]采用格子波爾茲曼方法（lattice Boltzmann method，LBM）模擬雷諾數(shù)在100～1000范圍的中性浮力顆粒在方形截面管道內(nèi)的遷移，通過計(jì)算顆粒邊界節(jié)點(diǎn)上動(dòng)量的變化得到顆粒受到的力和力矩，數(shù)值模擬得到的方形截面管道內(nèi)單個(gè)顆粒遷移的平衡位置與實(shí)驗(yàn)報(bào)道[1]的圓管內(nèi)顆粒遷移的平衡位置基本一致。研究發(fā)現(xiàn)，受雷諾數(shù)的影響，顆?？梢赃w移到拐角或邊界中心的若干個(gè)平衡位置。2009年，Gupta等[25]利用LBM模擬了矩形截面通道內(nèi)的中性浮力顆粒的慣性遷移，發(fā)現(xiàn)顆粒平衡位置受到雷諾數(shù)和通道橫截面高寬比的影響，并將模擬結(jié)果與已有文獻(xiàn)報(bào)道的實(shí)驗(yàn)結(jié)果[1]進(jìn)行對(duì)比分析。

2019年，Liu等[26]采用LBM-DEM（discrete element method，DEM）耦合對(duì)二維平面泊肅葉流下中性浮力顆粒的橫向遷移進(jìn)行數(shù)值模擬，研究發(fā)現(xiàn)，顆粒濃度對(duì)橫向遷移特性的影響主要取決于雷諾數(shù)Re，Re<20時(shí)，在顆粒體積分?jǐn)?shù)不超過5%時(shí)才能觀察到顆粒的遷移現(xiàn)象，Re>20時(shí)，在顆粒體積分?jǐn)?shù)高于20%時(shí)仍然能夠發(fā)生顆粒遷移，并提出了富集數(shù)以表征顆粒慣性遷移的程度。數(shù)值模擬得到的單個(gè)顆粒橫向遷移的無量綱平衡位置隨雷諾數(shù)的變化關(guān)系與實(shí)驗(yàn)結(jié)果[1]吻合良好。

2.4 界面追蹤方法

2011年，Nourbakhsh等[27]采用有限差分/界面追蹤方法模擬了三維可變形液滴在平面泊肅葉流中的運(yùn)動(dòng)，研究了毛細(xì)數(shù)（黏性力與界面張力的比值）、雷諾數(shù)、體積分?jǐn)?shù)對(duì)液滴最終穩(wěn)定位置的影響，數(shù)值模擬得到的通道截面上顆粒濃度分布與已有實(shí)驗(yàn)基本一致。2013年，Khalili等[28]利用界面追蹤方法數(shù)值模擬了泊肅葉流中二維浮力液滴的慣性遷移，研究了弗勞德數(shù)（浮力與慣性力的比值）、雷諾數(shù)、毛細(xì)數(shù)、液滴與流體的密度比對(duì)液滴最終平衡位置的影響，然而模擬結(jié)果未得到實(shí)驗(yàn)的證實(shí)。2013年，Bayareh等[29]利用有限差分/界面追蹤方法模擬了泊肅葉流中可變形非中性浮力液滴發(fā)生慣性遷移的平衡位置，研究發(fā)現(xiàn)，對(duì)于輕微浮力液滴，平衡位置在壁面和管道中心線之間。如果液滴滯后于流體，平衡位置接近中心線；如果液滴超前于流體，則平衡位置接近壁面。隨著雷諾數(shù)的增加，較輕液滴的平衡位置稍微接近于壁面，較重液滴的平衡位置則更接近中心線，但是，模擬預(yù)測結(jié)果仍有待實(shí)驗(yàn)驗(yàn)證。

2.5 其他方法

2005年，Yang等[20]采用任意拉格朗日?歐拉動(dòng)網(wǎng)格技術(shù)描述顆粒的運(yùn)動(dòng)，模擬了管內(nèi)單個(gè)中性浮力剛性顆粒的橫向遷移，得到的顆粒平衡位置與實(shí)驗(yàn)吻合良好。

2005年，Cho等[30]利用納維?斯托克斯方程與顆粒運(yùn)動(dòng)方程耦合的直接數(shù)值模擬方法研究二維通道內(nèi)顆粒負(fù)載流動(dòng)，證實(shí)了已有實(shí)驗(yàn)中觀察到的顆粒橫向遷移現(xiàn)象，研究了通道尺寸與顆粒粒徑之比對(duì)顆粒平衡位置的影響。

2015年，Pazouki等[31]在拉格朗日?拉格朗日方法框架下，利用光滑粒子流體動(dòng)力學(xué)方法（smoothed particle hydrodynamics,SPH）描述流場，利用牛頓?歐拉運(yùn)動(dòng)方程描述剛性顆粒的運(yùn)動(dòng)狀態(tài)，在較寬的雷諾數(shù)范圍內(nèi)借助已有的實(shí)驗(yàn)和數(shù)值模擬結(jié)果對(duì)模型進(jìn)行了驗(yàn)證，探討了顆粒旋轉(zhuǎn)、顆粒粒徑、顆粒間距、顆粒形狀（球形、橢球）、雷諾數(shù)對(duì)圓管泊肅葉流下中性浮力顆粒橫向遷移最終穩(wěn)定位置的影響。

2015年，Kim等[32]采用懲罰浸沒邊界方法（penalty immersed boundary method，PIBM）數(shù)值模擬了三維彈性膠囊在平面泊肅葉流下橫向遷移的平衡位置受顆粒的初始橫向位置、雷諾數(shù)、通道寬高比及流體黏度等的影響，然而缺少實(shí)驗(yàn)結(jié)果來驗(yàn)證數(shù)值模擬的準(zhǔn)確性。

3 圓管泊肅葉流下中性浮力顆粒橫向遷移的三維數(shù)值模擬

雖然一些學(xué)者對(duì)泊肅葉流下中性浮力顆粒的橫向遷移特性進(jìn)行了實(shí)驗(yàn)和數(shù)值模擬研究，但是，實(shí)驗(yàn)研究多關(guān)注顆粒遷移后最終穩(wěn)定的橫向位置，對(duì)顆粒遷移過程的研究存在明顯不足。受實(shí)驗(yàn)測試手段的限制，通過實(shí)驗(yàn)研究難以獲得顆粒動(dòng)力學(xué)行為和顆粒周圍流場的細(xì)節(jié)信息，因而也難以揭示顆粒橫向遷移的內(nèi)在機(jī)理。雖然多種數(shù)值模擬方法在顆粒橫向遷移研究中得到應(yīng)用，但是，數(shù)值模擬多針對(duì)二維平面泊肅葉流動(dòng)下顆粒的橫向遷移，對(duì)三維管道內(nèi)顆粒橫向遷移的研究相對(duì)較少。此外，針對(duì)高雷諾數(shù)時(shí)（Re>1000）圓管內(nèi)顆粒橫向遷移最終平衡位置的個(gè)數(shù)和位置的研究還存在不足，迫切需要對(duì)此開展三維圓管內(nèi)顆粒橫向遷移的全過程研究?，F(xiàn)將借助開源流體?顆粒系統(tǒng)數(shù)值模擬軟件CFDEMcoupling[33]，基于DF/FD方法模擬再現(xiàn)圓管泊肅葉流下單個(gè)顆粒橫向遷移行為，以期為今后通過數(shù)值模擬方法進(jìn)行顆粒橫向遷移的精細(xì)動(dòng)力學(xué)行為及其內(nèi)在機(jī)理研究提供基礎(chǔ)。

3.1 數(shù)學(xué)模型

針對(duì)三維水平圓管泊肅葉流下單個(gè)中性浮力顆粒的慣性遷移進(jìn)行數(shù)值模擬研究，所研究液固兩相系統(tǒng)的示意圖如圖5所示。L為管長，u為流體速度。

圖5 數(shù)值模擬中液固兩相系統(tǒng)示意圖Fig.5 Schematic diagram of the liquid-solid two-phase system in numerical simulation

采用DF/FD方法進(jìn)行數(shù)值模擬。首先，忽略顆粒的存在，整個(gè)區(qū)域都作為流體域計(jì)算；然后，計(jì)算流體對(duì)顆粒的作用力，將該力以體積力的形式包含到顆粒運(yùn)動(dòng)模型中；最后，為了滿足守恒方程，對(duì)速度場和壓力場進(jìn)行了修正[34]。

管內(nèi)流體的連續(xù)性方程和動(dòng)量方程為[34-36]

式中：u為速度矢量；ρ為流體密度；p為流體壓力；t為時(shí)間；μ為流體動(dòng)力黏度。

采用牛頓?歐拉方法對(duì)單個(gè)中性浮力顆粒運(yùn)動(dòng)進(jìn)行建模[34,36]，則有

式中：m，I分別為顆粒的質(zhì)量和轉(zhuǎn)動(dòng)慣量；v，ω分別為顆粒速度和角速度；Ff為顆粒受到的流體力；R為顆粒中心指向顆粒表面的矢量。

通過對(duì)顆粒的浸沒邊界受力進(jìn)行積分以求解流體力Ff[34,36]。

式中，Ωs為顆粒邊界。

在流固耦合計(jì)算過程中，首先通過動(dòng)量方程（式（2））計(jì)算得到臨時(shí)速度場，然后將顆粒的速度施加到臨時(shí)速度場中，得到新速度場；為了滿足連續(xù)性方程（式（1）），利用泊松方程修正新速度場，得到修正后的速度場；將修正后的速度場引入到動(dòng)量方程中，對(duì)壓力場進(jìn)行修正[34,36]。

3.2 數(shù)值模擬條件和方法

數(shù)值模擬時(shí)采用和Matas等[1]的實(shí)驗(yàn)相同的流體條件（ρ=1050 kg/m3，μ=1.52×10?3m2/s，平均流速U=0.0124 m/s）、顆粒參數(shù)（d=0.9 mm，密度與流體密度相同）及管徑（D=8 mm）。

為了減少計(jì)算量，一方面，對(duì)進(jìn)口和出口施加周期性邊界條件，即從出口流出的流體和顆粒，以相同的速度從進(jìn)口流入，采用更短的管長（L=100 mm）模擬管長較長的管內(nèi)顆粒的行為[37]；另一方面，網(wǎng)格劃分時(shí)，采用靜態(tài)網(wǎng)格和動(dòng)態(tài)網(wǎng)格加密相結(jié)合的方法。具體網(wǎng)格劃分方法：首先對(duì)整個(gè)計(jì)算區(qū)域劃分靜態(tài)網(wǎng)格（圖6）；計(jì)算過程中，對(duì)于流體與顆粒界面所在的靜態(tài)網(wǎng)格進(jìn)行動(dòng)態(tài)網(wǎng)格加密（圖7）。從而實(shí)現(xiàn)了對(duì)流固界面的精準(zhǔn)描述，同時(shí)避免了全部采用靜態(tài)網(wǎng)格時(shí)網(wǎng)格數(shù)目龐大的問題。

圖6 數(shù)值模擬采用的靜態(tài)網(wǎng)格Fig.6 Static mesh used in numerical simulation

圖7 顆粒與流體界面的靜態(tài)網(wǎng)格和動(dòng)態(tài)網(wǎng)格Fig. 7 Static and dynamic meshes at particle-fluid boundary used in numerical simulation

數(shù)值模擬中，采用PISO算法求解流場。鑒于求解流場的時(shí)間步長應(yīng)小于流體流經(jīng)一個(gè)網(wǎng)格的時(shí)間，即滿足庫朗（Courant）數(shù)小于1；求解顆粒運(yùn)動(dòng)的時(shí)間步長應(yīng)小于顆粒的松弛時(shí)間，結(jié)合流固耦合計(jì)算的特點(diǎn)，流場計(jì)算時(shí)間步長采用3×10?5s，顆粒計(jì)算時(shí)間步長采用3×10?6s。

3.3 泊肅葉流速度分布和顆粒橫向遷移過程數(shù)值模擬結(jié)果

數(shù)值模擬得到的遠(yuǎn)離顆粒位置（與顆粒的軸向間距大于10d）的管道橫截面上流體速度曲線與泊肅葉流下流體速度的解析解的對(duì)比關(guān)系如圖8所示，可見本文模擬得到的數(shù)值解和解析解吻合良好。u為當(dāng)?shù)亓黧w的速度，umax為圓管內(nèi)的最大速度。此外，從數(shù)值模擬結(jié)果可以看出，數(shù)值模擬和解析解得到的管道內(nèi)流體的速度分布幾乎完全重合，這表明建立的模型和計(jì)算方法能夠準(zhǔn)確描述泊肅葉流的流場特性。

圖8 泊肅葉流速度分布數(shù)值解與解析解的對(duì)比Fig.8 Comparison between numerical and analytical solutions of the velocity profile in Poiseuille flow

圖9給出了數(shù)值模擬得到的單個(gè)顆粒橫向遷移隨時(shí)間的變化關(guān)系，如果不考慮顆粒旋轉(zhuǎn)（不考慮式（4）），仍能發(fā)生顆粒的橫向遷移，無量綱初始位置r0/R=0.2，0.4，0.75時(shí)，顆粒最終穩(wěn)定位置對(duì)應(yīng)的r/R≈0.4，這與Matas等[1]的實(shí)驗(yàn)結(jié)果中顆粒橫向遷移最終穩(wěn)定的平衡位置r/R≈0.68仍有很大差異。在考慮顆粒旋轉(zhuǎn)（增加式（4））后，數(shù)值模擬得到的顆粒橫向遷移最終穩(wěn)定的平衡位置r/R≈0.66，數(shù)值模擬與實(shí)驗(yàn)的相對(duì)誤差小于3%。這表明在顆粒慣性遷移研究中應(yīng)當(dāng)考慮顆粒的旋轉(zhuǎn)。考慮旋轉(zhuǎn)時(shí)顆粒橫向遷移的平衡位置更接近壁面的原因是：顆粒運(yùn)動(dòng)對(duì)流體產(chǎn)生影響的距離對(duì)顆粒遷移的平衡位置起重要作用，該距離越大，顆粒越靠近壁面；旋轉(zhuǎn)顆粒運(yùn)動(dòng)比不旋轉(zhuǎn)顆粒運(yùn)動(dòng)對(duì)流體產(chǎn)生影響的距離更大，使得顆粒遷移的平衡位置更接近壁面。圖9的結(jié)果還表明，本文的模型和數(shù)值模擬方法能夠準(zhǔn)確預(yù)測圓管泊肅葉流下中性浮力顆粒的橫向遷移。

圖9 數(shù)值模擬得到的顆粒橫向遷移隨時(shí)間的變化關(guān)系Fig.9 Lateral migration of the particle as a function of time obtained from numerical simulation

4 結(jié)束語

隨著微流控技術(shù)的發(fā)展，管道內(nèi)顆粒橫向遷移現(xiàn)象成為研究熱點(diǎn)，一些學(xué)者針對(duì)泊肅葉流下中性浮力顆粒的橫向遷移開展了實(shí)驗(yàn)與數(shù)值模擬研究。然而，已有實(shí)驗(yàn)研究主要集中在對(duì)不同條件下顆粒橫向遷移最終穩(wěn)定平衡位置的研究，受限于微尺度條件下的測量手段，單純依靠實(shí)驗(yàn)研究難以展示顆粒橫向遷移的全過程及其內(nèi)在機(jī)理。數(shù)值模擬方法易于改變參數(shù)，便于獲取流場和顆粒動(dòng)力學(xué)特性的詳細(xì)信息，在揭示兩相流現(xiàn)象及機(jī)理上具有實(shí)驗(yàn)研究難以比擬的優(yōu)點(diǎn)。目前，有限元方法、虛擬區(qū)域方法、格子波爾茲曼方法、界面追蹤方法等先后被用于研究泊肅葉流下顆粒橫向遷移行為及機(jī)理。然而，由于顆粒橫向遷移中流場和顆粒參數(shù)的復(fù)雜性，迄今對(duì)顆粒橫向遷移行為的認(rèn)識(shí)仍不夠全面，甚至存在爭議之處，對(duì)其內(nèi)在機(jī)理的掌握還不夠深入。鑒于上述原因，本文借助開源流體?顆粒系統(tǒng)數(shù)值模擬軟件CFDEMcoupling，基于直接力/虛擬區(qū)域方法，開展泊肅葉流下單個(gè)中性浮力顆粒橫向遷移的數(shù)值模擬，獲得了準(zhǔn)確的泊肅葉流場信息，展現(xiàn)了顆粒橫向遷移的過程，為探究顆粒橫向遷移動(dòng)力學(xué)行為及機(jī)理提供了重要基礎(chǔ)。在今后的研究中，應(yīng)基于直接力/虛擬區(qū)域方法繼續(xù)開展研究，以全面識(shí)別顆粒橫向遷移的動(dòng)力學(xué)行為，并揭示其內(nèi)在機(jī)理，以期為顆粒橫向遷移的實(shí)際應(yīng)用提供科學(xué)依據(jù)和方法指導(dǎo)。