葉 挺,凡鳳仙
(上海理工大學 能源與動力工程學院,上海 200093)
管道內顆粒負載流動是能源、化工及生命科學等領域的研究熱點之一。顆粒橫向遷移是管道內顆粒負載流動中的一種有趣現象,它是指隨機散布在泊肅葉流中的中性浮力顆粒(密度與流體密度相同的顆粒),經過一定時間后,遷移到距管中心一定距離橫向位置上流動的現象[1]。顆粒橫向遷移現象有望在顆粒的分離、篩選、混合、反應中得到應用[2-3],富有研究價值和應用潛力。顆粒橫向遷移的早期研究主要采用實驗方法,且多關注雷諾數、顆粒尺寸等對顆粒遷移的最終平衡位置的影響。隨著計算機和數值模擬技術的發展,數值模擬逐漸成為研究顆粒負載流動的重要手段,通過數值模擬可以方便地對參數進行調整,且能夠獲取豐富的流場和顆粒動態過程的細節信息,便于揭示相關現象的規律及內在機理。近年來,多種數值模擬方法在顆粒橫向遷移研究中得到應用,并取得了一些有價值的研究成果,但由于顆粒橫向遷移所涉及的液固兩相流場的復雜性,目前對顆粒橫向遷移行為的認識仍存在爭議,對橫向遷移內在機理研究還不夠深入。本文首先針對泊肅葉流下中性浮力顆粒橫向遷移的相關實驗和數值模擬研究進行歸納和評述,然后基于直接力/虛擬區域(direct-forcing fictitious domain,DF/FD)方法模擬再現單個顆粒橫向遷移的過程,在此基礎上指出今后進行顆粒橫向遷移研究的思路和方向。
1961年,Segré和Silberberg[4]首次在實驗中觀察到圓管內低雷諾數(7 圖1 橫向遷移后顆粒濃度沿徑向的分布Fig. 1 Particle concentration as a function of radial position after lateral migration 隨后,一些學者迅速開展了不同條件下圓管內顆粒橫向遷移平衡位置的探討。1962年,Oliver[5]通過實驗發現顆粒的旋轉對于顆粒遷移的最終穩定位置有顯著影響,無旋轉時顆粒最終穩定位置更接近管中心軸線,并認為向內的壁面效應與向外的馬格努斯效應使得顆粒最終穩定在特定的徑向位置。1965年,Jeffrey等[6]在圓管直徑D=32.5 mm,顆粒直徑d=1.5,2.0,2.9 mm的條件下,實驗觀察到中性浮力顆粒最終遷移到距管中心軸線約為2R/3的環形區域。1966年,Karnis等[7]通過實驗發現中性浮力剛性顆粒橫向遷移的平衡位置與顆粒形狀、粒徑和管徑比有關,中性浮力可變形顆粒總是遷移到管中心軸線位置。 然而,由于受實驗技術的限制,難以得到更為詳細的顆粒動力學信息和流場信息,相關實驗研究一度停滯。直到進入21世紀,微流控技術興起,泊肅葉流下中性浮力顆粒的橫向遷移重新引起了學者的重視。2004年,Eloot等[8]實驗研究了顆粒直徑d=10 μm的近中性浮力球形顆粒在直徑D=220,530 μm的毛細管內的運動,結果表明,在較高的雷諾數、較長和較細的毛細管、較稀的顆粒濃度的條件下,顆粒橫向遷移的平衡位置更接近于管壁。2004年,Matas等[1]實驗研究了雷諾數Re在67~1700范圍以及達到2 000時,圓管內稀相泊肅葉流下中性浮力球形顆粒(D/d=8~42)的橫向遷移,發現Re較小時,顆粒遷移到距管中心軸線一定距離的環形區域,隨著Re的增大,環形區域的位置向管壁移動,當Re繼續增大到一定值時,顆粒的遷移使得顆粒分布在2個環形區域,即除了靠近管壁的外環外,還存在一個靠近管中心的內環。圖2給出了Matas等[1]實驗得到的Re=1000時顆粒概率密度P沿徑向的分布情況,可見內環出現在r/R≈0.5處,外環出現在r/R≈0.9處。2017年,Morita等[9]實驗研究了雷諾數Re在100~1100范圍、顆粒粒徑與管徑比約為0.1時,泊肅葉流下中性浮力顆粒橫向遷移的平衡位置,研究發現,Re較大時,顆粒分布首先集中在2個環形(即內環和外環)區域,隨著顆粒向下游運動,內環附近顆粒減少,而外環附近顆粒增加,若管道足夠長,顆粒將最終遷移到外環位置。Morita等[9]得到了較大雷諾數下顆粒遷移特性的演變過程示意圖,如圖3所示,即隨機分布的顆粒隨流體進入管內(x=0),首先向內遷移到達一定的位置(x=L2);接著轉而向外遷移,引起顆粒在內環和外環2個區域集中(x=L1);顆粒繼續向外遷移,最終集中在一個環形區域(x=L0)。在圖3中,顆粒向內遷移發生在流體從初始狀態充分發展為泊肅葉流的管段,這表明顆粒起初向內遷移與流動的發展過程有關,顆粒轉而向外遷移則是由于泊肅葉流中橫向速度梯度引起的升力和顆粒的旋轉引起的升力所致。 圖2 橫向遷移后顆粒概率密度沿徑向的分布Fig.2 Particle probability density as a function of radial position after lateral migration 圖3 較大雷諾數下顆粒遷移特性的演變過程示意圖Fig.3 Schematic diagram of the evolution of particle migration at relatively high Reynolds number 由于壁面效應的影響,非中心對稱的矩形和方形截面通道內顆粒橫向遷移的特性將與中心對稱的圓管不同,一些學者針對矩形與方形截面通道內顆粒的橫向遷移開展了一系列研究。 在矩形截面通道泊肅葉流下顆粒橫向遷移的研究方面,2005年,Staben等[10]實驗研究了雷諾數Re(以矩形窄邊長H作為特征尺寸)在0.005~0.05范圍時顆粒粒徑對不同管口形狀的矩形截面微通道內顆粒橫向遷移的影響,結果表明,小粒徑顆粒(0.07≤d/H≤0.10)在直管口通道內向通道中心遷移,在帶倒角管口通道內分布更為均勻;中等尺度粒徑顆粒(0.46≤d/H≤0.54)在直管口和帶倒角管口通道內均呈幾乎均勻的分布;少量大粒徑顆粒(0.79≤d/H≤0.95)的分布略微向通道的一個側壁傾斜。2008年,Bhagat等[11]在Re>1的條件下實驗研究了矩形截面微通道內顆粒的橫向遷移特性,結果表明,顆粒遷移過程主要取決于通道截面較短邊的幾何尺度,顆粒最終分布在較大的壁面附近,d/H=0.07情況下,Re≥ 10時才能發生橫向遷移,且隨著d/H的增加,顆粒橫向遷移所需要的臨界雷諾數降低。由于文獻[10-11]的實驗中采用的雷諾數范圍有很大差異,使得實驗結果不同,雷諾數對矩形微通道內顆粒橫向遷移的影響仍需要進一步系統的研究。 在方形截面通道泊肅葉流下顆粒橫向遷移的研究方面,2008年,Kim等[12]實驗研究了通道雷諾數Re在0.06~58.65范圍內中性浮力顆粒(通道水力直徑與粒徑之比λ≈14)的橫向遷移,結果表明,在很低的雷諾數下,顆粒的橫向遷移仍會發生,存在一個在20~30之間的臨界雷諾數Rec,當Re 圖4 典型雷諾數下橫向遷移后的顆粒分布Fig.4 Particle distributions after lateral migration at typical Reynolds numbers 為了理解顆粒橫向遷移的機理,1979年,Aokl等[15]實驗研究了基于顆粒自由沉降速度的雷諾數ReF(ReF=2r'uF/ν,r'為顆粒半徑,uF為顆粒在靜止流體中的自由沉降速度,ν為流體運動黏度)在0.043~10范圍時非中性浮力顆粒在圓管內的橫向遷移,結果表明,ReF>1時,若顆粒運動滯后于流體,則顆粒向管中心軸線遷移,若顆粒運動超前于流體,則顆粒向管壁遷移,這一現象受控于Magnus效應;ReF<1時,靠近管壁的顆粒向內遷移,靠近管中心的顆粒向外遷移,最終顆粒到達管壁和管中心軸線之間的一個平衡位置,這一現象是由壁面效應和慣性作用共同引起的。為了探討顆粒橫向遷移在微流控領域的應用,2011年,文獻[16]基于螺旋微通道內慣性升力和Dean曳力共同作用下,不同粒徑顆粒發生橫向遷移的最終平衡位置不同的原理,設計出利用螺旋形通道富集和分離顆粒的微流控裝置。2012年,Martel等[17]實驗研究了高寬比在0.125~1范圍的螺旋形微通道內顆粒富集動力學行為,給出了實現高質量富集所需要的通道長度和流速范圍。這些研究可以為理解顆粒橫向遷移的機理以及開發基于橫向遷移原理的微流控裝置提供參考。 與實驗研究相比,管道內中性浮力顆粒橫向遷移的數值模擬研究開展得較晚,但由于數值模擬的迅速發展以及計算機能力的快速提升,數值模擬已成為研究顆粒橫向遷移行為及其內在機理的重要手段。目前,研究顆粒橫向遷移的數值模擬方法主要有有限元方法、虛擬區域方法、格子波爾茲曼方法及界面追蹤方法等。 1994年,Feng等[18]對泊肅葉流下中性浮力圓形顆粒的運動進行二維有限元模擬,獲得了與文獻[4]的實驗類似的顆粒橫向遷移效應,并將顆粒橫向遷移的驅動力歸因于潤滑產生的壁面排斥力、剪切滑移產生的慣性升力、顆粒旋轉產生的升力和速度曲線的曲率產生的升力。 2009年,文獻[19]利用有限元方法求解穩態納維?斯托克斯方程,通過表面積分的方法獲得顆粒受到的力和力矩,對狹窄矩形截面微通道內顆粒橫向遷移進行數值模擬,并進行了相應的實驗,數值模擬得到的顆粒橫向遷移的平衡位置與實驗結果吻合良好。在此基礎上,模擬了顆粒橫向遷移過程中的升力,分析了升力的標度關系。 2005年,Yang等[20]采用虛擬區域方法,通過固定網格分布式拉格朗日乘子描述顆粒的運動,模擬了管內單個中性浮力剛性顆粒的橫向遷移,得到的顆粒平衡位置與實驗吻合良好。2008年,Pan等[21]采用基于拉格朗日乘子的虛擬區域方法模擬三維泊肅葉流中的中性浮力橢圓形顆粒的運動,研究了顆粒的旋轉和取向特性,發現隨著雷諾數和顆粒形狀的變化,顆粒呈現截然不同的運動狀態,顆粒橫向遷移的平衡位置與文獻報道的圓盤形顆粒的實驗結果[7]類似。上述方法又稱為DLM/FD(distributed-Lagrange-multiplier fictitious domain)方法,該方法需要顯式計算拉格朗日乘子進而計算顆粒的運動參數,因而計算效率不高。 為了提高虛擬區域方法的計算效率,直接力/虛擬區域(direct-forcing fictitious domain,DF/FD)方法被提出,該方法引入體積力以計算顆粒的運動參數,從而避免了對拉格朗日乘子的依賴。2008年,Shao等[22]在流動方向上施加周期性邊界條件,利用DF/FD方法模擬了圓管泊肅葉流下球形中性浮力顆粒的運動過程以及顆粒運動對流場的影響,發現在雷諾數高于臨界值時顆粒橫向遷移到一定的平衡位置,這與文獻[1]的實驗結果在定性上一致。2009年,Sun等[23]采用DF/FD方法模擬了二維通道內的圓形中性浮力顆粒在恒定壓力和波動壓力驅動的流動中,橫向遷移過程以及顆粒附近的流場狀況,但缺少和實驗結果的對比驗證。 2006年,Chun等[24]采用格子波爾茲曼方法(lattice Boltzmann method,LBM)模擬雷諾數在100~1000范圍的中性浮力顆粒在方形截面管道內的遷移,通過計算顆粒邊界節點上動量的變化得到顆粒受到的力和力矩,數值模擬得到的方形截面管道內單個顆粒遷移的平衡位置與實驗報道[1]的圓管內顆粒遷移的平衡位置基本一致。研究發現,受雷諾數的影響,顆??梢赃w移到拐角或邊界中心的若干個平衡位置。2009年,Gupta等[25]利用LBM模擬了矩形截面通道內的中性浮力顆粒的慣性遷移,發現顆粒平衡位置受到雷諾數和通道橫截面高寬比的影響,并將模擬結果與已有文獻報道的實驗結果[1]進行對比分析。 2019年,Liu等[26]采用LBM-DEM(discrete element method,DEM)耦合對二維平面泊肅葉流下中性浮力顆粒的橫向遷移進行數值模擬,研究發現,顆粒濃度對橫向遷移特性的影響主要取決于雷諾數Re,Re<20時,在顆粒體積分數不超過5%時才能觀察到顆粒的遷移現象,Re>20時,在顆粒體積分數高于20%時仍然能夠發生顆粒遷移,并提出了富集數以表征顆粒慣性遷移的程度。數值模擬得到的單個顆粒橫向遷移的無量綱平衡位置隨雷諾數的變化關系與實驗結果[1]吻合良好。 2011年,Nourbakhsh等[27]采用有限差分/界面追蹤方法模擬了三維可變形液滴在平面泊肅葉流中的運動,研究了毛細數(黏性力與界面張力的比值)、雷諾數、體積分數對液滴最終穩定位置的影響,數值模擬得到的通道截面上顆粒濃度分布與已有實驗基本一致。2013年,Khalili等[28]利用界面追蹤方法數值模擬了泊肅葉流中二維浮力液滴的慣性遷移,研究了弗勞德數(浮力與慣性力的比值)、雷諾數、毛細數、液滴與流體的密度比對液滴最終平衡位置的影響,然而模擬結果未得到實驗的證實。2013年,Bayareh等[29]利用有限差分/界面追蹤方法模擬了泊肅葉流中可變形非中性浮力液滴發生慣性遷移的平衡位置,研究發現,對于輕微浮力液滴,平衡位置在壁面和管道中心線之間。如果液滴滯后于流體,平衡位置接近中心線;如果液滴超前于流體,則平衡位置接近壁面。隨著雷諾數的增加,較輕液滴的平衡位置稍微接近于壁面,較重液滴的平衡位置則更接近中心線,但是,模擬預測結果仍有待實驗驗證。 2005年,Yang等[20]采用任意拉格朗日?歐拉動網格技術描述顆粒的運動,模擬了管內單個中性浮力剛性顆粒的橫向遷移,得到的顆粒平衡位置與實驗吻合良好。 2005年,Cho等[30]利用納維?斯托克斯方程與顆粒運動方程耦合的直接數值模擬方法研究二維通道內顆粒負載流動,證實了已有實驗中觀察到的顆粒橫向遷移現象,研究了通道尺寸與顆粒粒徑之比對顆粒平衡位置的影響。 2015年,Pazouki等[31]在拉格朗日?拉格朗日方法框架下,利用光滑粒子流體動力學方法(smoothed particle hydrodynamics,SPH)描述流場,利用牛頓?歐拉運動方程描述剛性顆粒的運動狀態,在較寬的雷諾數范圍內借助已有的實驗和數值模擬結果對模型進行了驗證,探討了顆粒旋轉、顆粒粒徑、顆粒間距、顆粒形狀(球形、橢球)、雷諾數對圓管泊肅葉流下中性浮力顆粒橫向遷移最終穩定位置的影響。 2015年,Kim等[32]采用懲罰浸沒邊界方法(penalty immersed boundary method,PIBM)數值模擬了三維彈性膠囊在平面泊肅葉流下橫向遷移的平衡位置受顆粒的初始橫向位置、雷諾數、通道寬高比及流體黏度等的影響,然而缺少實驗結果來驗證數值模擬的準確性。 雖然一些學者對泊肅葉流下中性浮力顆粒的橫向遷移特性進行了實驗和數值模擬研究,但是,實驗研究多關注顆粒遷移后最終穩定的橫向位置,對顆粒遷移過程的研究存在明顯不足。受實驗測試手段的限制,通過實驗研究難以獲得顆粒動力學行為和顆粒周圍流場的細節信息,因而也難以揭示顆粒橫向遷移的內在機理。雖然多種數值模擬方法在顆粒橫向遷移研究中得到應用,但是,數值模擬多針對二維平面泊肅葉流動下顆粒的橫向遷移,對三維管道內顆粒橫向遷移的研究相對較少。此外,針對高雷諾數時(Re>1000)圓管內顆粒橫向遷移最終平衡位置的個數和位置的研究還存在不足,迫切需要對此開展三維圓管內顆粒橫向遷移的全過程研究?,F將借助開源流體?顆粒系統數值模擬軟件CFDEMcoupling[33],基于DF/FD方法模擬再現圓管泊肅葉流下單個顆粒橫向遷移行為,以期為今后通過數值模擬方法進行顆粒橫向遷移的精細動力學行為及其內在機理研究提供基礎。 針對三維水平圓管泊肅葉流下單個中性浮力顆粒的慣性遷移進行數值模擬研究,所研究液固兩相系統的示意圖如圖5所示。L為管長,u為流體速度。 圖5 數值模擬中液固兩相系統示意圖Fig.5 Schematic diagram of the liquid-solid two-phase system in numerical simulation 采用DF/FD方法進行數值模擬。首先,忽略顆粒的存在,整個區域都作為流體域計算;然后,計算流體對顆粒的作用力,將該力以體積力的形式包含到顆粒運動模型中;最后,為了滿足守恒方程,對速度場和壓力場進行了修正[34]。 管內流體的連續性方程和動量方程為[34-36] 式中:u為速度矢量;ρ為流體密度;p為流體壓力;t為時間;μ為流體動力黏度。 采用牛頓?歐拉方法對單個中性浮力顆粒運動進行建模[34,36],則有 式中:m,I分別為顆粒的質量和轉動慣量;v,ω分別為顆粒速度和角速度;Ff為顆粒受到的流體力;R為顆粒中心指向顆粒表面的矢量。 通過對顆粒的浸沒邊界受力進行積分以求解流體力Ff[34,36]。 式中,Ωs為顆粒邊界。 在流固耦合計算過程中,首先通過動量方程(式(2))計算得到臨時速度場,然后將顆粒的速度施加到臨時速度場中,得到新速度場;為了滿足連續性方程(式(1)),利用泊松方程修正新速度場,得到修正后的速度場;將修正后的速度場引入到動量方程中,對壓力場進行修正[34,36]。 數值模擬時采用和Matas等[1]的實驗相同的流體條件(ρ=1050 kg/m3,μ=1.52×10?3m2/s,平均流速U=0.0124 m/s)、顆粒參數(d=0.9 mm,密度與流體密度相同)及管徑(D=8 mm)。 為了減少計算量,一方面,對進口和出口施加周期性邊界條件,即從出口流出的流體和顆粒,以相同的速度從進口流入,采用更短的管長(L=100 mm)模擬管長較長的管內顆粒的行為[37];另一方面,網格劃分時,采用靜態網格和動態網格加密相結合的方法。具體網格劃分方法:首先對整個計算區域劃分靜態網格(圖6);計算過程中,對于流體與顆粒界面所在的靜態網格進行動態網格加密(圖7)。從而實現了對流固界面的精準描述,同時避免了全部采用靜態網格時網格數目龐大的問題。 圖6 數值模擬采用的靜態網格Fig.6 Static mesh used in numerical simulation 圖7 顆粒與流體界面的靜態網格和動態網格Fig. 7 Static and dynamic meshes at particle-fluid boundary used in numerical simulation 數值模擬中,采用PISO算法求解流場。鑒于求解流場的時間步長應小于流體流經一個網格的時間,即滿足庫朗(Courant)數小于1;求解顆粒運動的時間步長應小于顆粒的松弛時間,結合流固耦合計算的特點,流場計算時間步長采用3×10?5s,顆粒計算時間步長采用3×10?6s。 數值模擬得到的遠離顆粒位置(與顆粒的軸向間距大于10d)的管道橫截面上流體速度曲線與泊肅葉流下流體速度的解析解的對比關系如圖8所示,可見本文模擬得到的數值解和解析解吻合良好。u為當地流體的速度,umax為圓管內的最大速度。此外,從數值模擬結果可以看出,數值模擬和解析解得到的管道內流體的速度分布幾乎完全重合,這表明建立的模型和計算方法能夠準確描述泊肅葉流的流場特性。 圖8 泊肅葉流速度分布數值解與解析解的對比Fig.8 Comparison between numerical and analytical solutions of the velocity profile in Poiseuille flow 圖9給出了數值模擬得到的單個顆粒橫向遷移隨時間的變化關系,如果不考慮顆粒旋轉(不考慮式(4)),仍能發生顆粒的橫向遷移,無量綱初始位置r0/R=0.2,0.4,0.75時,顆粒最終穩定位置對應的r/R≈0.4,這與Matas等[1]的實驗結果中顆粒橫向遷移最終穩定的平衡位置r/R≈0.68仍有很大差異。在考慮顆粒旋轉(增加式(4))后,數值模擬得到的顆粒橫向遷移最終穩定的平衡位置r/R≈0.66,數值模擬與實驗的相對誤差小于3%。這表明在顆粒慣性遷移研究中應當考慮顆粒的旋轉。考慮旋轉時顆粒橫向遷移的平衡位置更接近壁面的原因是:顆粒運動對流體產生影響的距離對顆粒遷移的平衡位置起重要作用,該距離越大,顆粒越靠近壁面;旋轉顆粒運動比不旋轉顆粒運動對流體產生影響的距離更大,使得顆粒遷移的平衡位置更接近壁面。圖9的結果還表明,本文的模型和數值模擬方法能夠準確預測圓管泊肅葉流下中性浮力顆粒的橫向遷移。 圖9 數值模擬得到的顆粒橫向遷移隨時間的變化關系Fig.9 Lateral migration of the particle as a function of time obtained from numerical simulation 隨著微流控技術的發展,管道內顆粒橫向遷移現象成為研究熱點,一些學者針對泊肅葉流下中性浮力顆粒的橫向遷移開展了實驗與數值模擬研究。然而,已有實驗研究主要集中在對不同條件下顆粒橫向遷移最終穩定平衡位置的研究,受限于微尺度條件下的測量手段,單純依靠實驗研究難以展示顆粒橫向遷移的全過程及其內在機理。數值模擬方法易于改變參數,便于獲取流場和顆粒動力學特性的詳細信息,在揭示兩相流現象及機理上具有實驗研究難以比擬的優點。目前,有限元方法、虛擬區域方法、格子波爾茲曼方法、界面追蹤方法等先后被用于研究泊肅葉流下顆粒橫向遷移行為及機理。然而,由于顆粒橫向遷移中流場和顆粒參數的復雜性,迄今對顆粒橫向遷移行為的認識仍不夠全面,甚至存在爭議之處,對其內在機理的掌握還不夠深入。鑒于上述原因,本文借助開源流體?顆粒系統數值模擬軟件CFDEMcoupling,基于直接力/虛擬區域方法,開展泊肅葉流下單個中性浮力顆粒橫向遷移的數值模擬,獲得了準確的泊肅葉流場信息,展現了顆粒橫向遷移的過程,為探究顆粒橫向遷移動力學行為及機理提供了重要基礎。在今后的研究中,應基于直接力/虛擬區域方法繼續開展研究,以全面識別顆粒橫向遷移的動力學行為,并揭示其內在機理,以期為顆粒橫向遷移的實際應用提供科學依據和方法指導。


1.2 矩形與方形通道內的實驗研究

1.3 其他相關研究
2 泊肅葉流下中性浮力顆粒橫向遷移的數值模擬研究
2.1 有限元方法
2.2 虛擬區域方法
2.3 格子玻爾茲曼方法
2.4 界面追蹤方法
2.5 其他方法
3 圓管泊肅葉流下中性浮力顆粒橫向遷移的三維數值模擬
3.1 數學模型




3.2 數值模擬條件和方法


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


4 結束語