方管中中性懸浮顆粒慣性遷移現象的數值研究*

姚同林，余釗圣，邵雪明

(浙江大學力學系，浙江 杭州 310027)

0 引 言

顆粒懸浮流在自然界和工農業生產中廣泛存在，如流化床、管道輸運、細胞分離，以及更為常見的河道中泥沙沉降和輸運等。盡管學術界對這種流動已有了大量的實驗、數值方面的研究，但由于顆粒與流體之間相互作用的復雜性，至今對其流動特性以及作用機理仍缺乏足夠的認識。

其中，Poiseuille流動中的粒子遷移問題是一個研究熱點。該問題起源于 Segré and Silberberg[1－2]在1961年觀測到的中性懸浮圓形顆粒在圓管中遷移到一個平衡位置，該平衡位置大約距離圓管中心軸有0.6倍的管半徑的距離。之后，更多的學者(如Jeffrey[3]，Karnis[4]，Matas[5－6]等)通過更多深入的實驗研究確認了Segré—Silberberg效應。最近，Matas等在 Re高于2400(基于管直徑和平均速度)的實驗中，發現了更靠近中心軸的平衡位置，稱之為內環。在粒子遷移問題的數值模擬中，J.Feng等[7]研究了單個圓形粒子在二維 Poiseuille 流動中的運動。Pan ＆ Glowinski[8－9]在平板和圓管Poiseuille流動下模擬了中性懸浮粒子的運動。Yu等[10]在垂直管 Poiseuille流動中模擬了中性懸浮和非中性懸浮球形顆粒的徑向遷移速度、角速度及軸向速度。Shao等[11]在較大Re下，對圓管中球形顆粒的運動進行了模擬，驗證了內部平衡點的存在。Chun＆Ladd[12]研究了球形顆粒在方管中的遷移現象，模擬出了顆粒的兩種平衡位置，并提出將兩個粒子用彈簧捆綁在一起形成啞鈴結構，新的平衡位置則會偏向方管中心。

綜上所述，文獻中對于圓管顆粒問題研究比較多，但對于方管中顆粒的慣性遷移問題模擬工作比較少，更缺乏對形成的兩種不同平衡位置的研究。

本研究采用虛擬區域方法，研究中性懸浮球形顆粒在上下前后邊界均為固定壁面的方管Poiseuille流中的慣性遷移問題。本研究主要關注雷諾數為100～1500范圍內，顆粒粒徑對慣性遷移現象的影響。本研究算法中對湍流和顆粒的模擬均使用直接數值模擬，沒有引入模型，所以稱為完全直接數值模擬。所采用的虛擬區域方法基于非貼體網格，在計算過程中無需移動和重新劃分網格，簡單高效且精度高。但對于網格大小以及時間步長有比較高的要求，而且計算量大。

1 數值方法

1.1 物理模型

一個顆粒在方管中的簡單示意圖如圖1所示，其中z表示方管的流向，球形顆粒的半徑為a;方管邊長為2H;方管長度為L。在數值模擬中，顆粒與流體密度比ρr=1，即在方管內球形粒子處于中性懸浮狀態。為減少計算區域，流向采用周期性邊界條件，從右邊界流出的顆粒與流體，再次以相同條件從左邊界進入計算區域，如此循環，即可代表無限長的管道;x，y方向均采用壁面無滑移邊界條件。

圖1 方管中一個顆粒的示意圖

本研究中將方管邊長的一半(H)和流向中心線速度(um)作為特征長度和特征速度，以此定義的雷諾數為Rem=umH/ν。由于流向采用周期性邊界條件，需要額外的壓力梯度▽p=－dp/dz克服壁面摩擦力，以維持流動。方管進口速度剖面為方管Poiseille流動速度剖面:

對處于平衡態的方管中的流體，無量綱后的額外壓力梯度為:

1.2 虛擬區域法

本研究采用的模擬方法是由Glowinski等人[13]提出、余釗圣[14]推廣改進的基于非貼體網格的直接力/虛擬區域法(DF/FD)。該方法的核心思想是假設固體顆粒內部充滿流體，通過在這部分流體上引入虛擬體積力，使其符合剛體顆粒的運動。基于上述假設，虛擬區域法很好地處理了流體跟固體顆粒之間的界面問題以及兩者之間的相互作用問題，在槽流和管流中得到了廣泛應用，其有效性和正確性得到了充分的證明。

基于虛擬區域法的假定，本研究在顆粒內部引入虛擬體積力λ，假設顆粒的密度、體積、轉動慣量、速度和角速度分別為 ρs，Vp，J，U 和 ωp，采用下述特征量來無量綱化控制方程:特征長度Lc，特征速度Uc，特征時間Lc/Uc，特征密度 ρf，特征壓力，特征擬體力。

無量綱化之后的控制方程為:

式中:Ω—整個計算區域;P(t)—顆粒占據的區域;Ωf—流體區域，Ωf=ΩP(t);r—以顆粒中心為原點的位置矢徑。

無量綱化控制參數如下:

密度比:ρr= ρs/ρf;

雷諾數:Re= ρfUcLc/μ;

通過時間分裂步格式離散上述控制方程，將原來流固耦合的問題式(3～7)分解成流體子問題和顆粒子問題。其中流體子問題是一個標準N－S方程的求解問題，利用基于半交錯風格的有限差分法和投影格式來求解，對于壓力泊松方程采用基于快速傅里葉變換技巧快速求解。對于空間離散，全部采用基于半交錯網格的二階精度有限差分格式。顆粒問題的求解與DLM/FD方法不同，不需要迭代求解。

另外，本研究利用區域分解和MPI實現算法的并行化。在并行算法中，速度方程和壓力方程均采用多重網格法(MG)進行迭代求解，有效保證了算法求解的速度與精度。

2 顆粒的慣性遷移分析

本節主要研究Re在100～1500之間，顆粒粒徑對顆粒在方管中的慣性遷移平衡位置的影響，分別計算了管長為2H條件下，顆粒尺寸a/H=0.15、a/H=0.1和a/H=0.075這3種工況，劃分的網格使得顆粒直徑上有10個節點，單個顆粒算例的具體參數如表1所示。每種工況下，筆者選取約20種具有代表性的顆粒的初始位置，具體布置如圖(2～5)所示。

表1 單個顆粒算例的具體參數

以a/H=0.15為例，圖(2～5)顯示了顆粒遷移軌跡，方管中顆粒的平衡位置有兩種:一種位于對角線上靠近角落;另一種位于邊線中部，在整個方管中存在8個平衡點。兩種平衡位置與壁面之間的距離有著明顯的區別，在Re數(如圖2、圖3所示)較小時，二者與壁面間距離相差不多，隨著Re數的增大，邊線中部的平衡位置明顯遠離壁面，而Chun＆Ladd采用格子玻爾茲曼方法模擬的結果顯示，兩個平衡位置總是距壁面0.2H。

圖2 Re=200下顆粒遷移軌跡在流向速度等值上的投影

圖3 Re=500下顆粒遷移軌跡在流向速度等值上的投影

圖4 Re=1000下顆粒遷移軌跡在流向速度等值上的投影

圖5 Re=1500下顆粒遷移軌跡在流向速度等值上的投影

a/H=0.15條件下，方管、圓管和槽道遷移平衡位置的比較如圖6所示。方管和圓管中的顆粒平衡位置相近，與槽道則明顯不同。顆粒的遷移，如前所述，源于固壁施加給顆粒的力與速度梯度產生的力之間的不平衡，而方管與圓管在兩個方向上為固壁邊界，而槽道只有法向為固壁邊界，這可以被認為是差異產生的主要原因。方管中對角線上的平衡位置，隨著Re數從100增加到1500，逐漸向管壁靠近，直至顆粒中心距管壁約為0.2H處，這與實驗中結果相似。方管邊線中間處的平衡位置，隨著Re數的增加，先靠近壁面，在Re數增加到800附近后，反而出現向方管中心遷移的現象，形成類似于圓管的內部平衡點。

圖6 a/H=0.15條件下，方管、圓管和槽道中顆粒遷移平衡位置的比較

兩種不同粒徑的顆粒都具有兩種平衡位置，且兩種位置都能夠穩定存在。隨著Re數的增加，3種大小顆粒的平衡位置改變的趨勢相同，如圖7所示。其區別主要體現在:相同Re數下，相對于大顆粒，小顆粒的平衡位置要更靠近于壁面。高Re數下，邊線中部的平衡位置相差較大，大顆粒更易形成內部平衡點。

圖7 a/H=0.1、a/H=0.15 和 a/H=0.075 顆粒的平衡位置

3 結束語

本研究采用并行虛擬區域方法對方管中中性懸浮顆粒慣性遷移現象進行了完全直接數值模擬，研究了雷諾數為100～1500范圍內，顆粒粒徑對慣性遷移現象的影響，豐富和完善了顆粒遷移運動的機理研究。研究結果表明:顆粒的遷移平衡位置分為方管對角線和邊線中間兩種。在周期性管長(L=2H)條件下，其平衡位置與粒徑大小密切相關。穩定存在的對角線和邊線中間兩個平衡位置，隨著雷諾數的增加，對角線上的平衡位置愈靠近角落，邊線中間的平衡位置則先靠近壁面，再遠離壁面，向管芯靠近;相同Re條件下，小顆粒的平衡位置比大顆粒更靠近壁面，尤其邊線中部的平衡位置差別較大。高Re條件下，大顆粒更易形成類似管流中的內部平衡點。

[1]SEGRE G，SILBERBERG A.Behaviour of macroscopic rigid spheres in poiseuille flow part 2.experimental results and interpretation[J].Journal of Fluid Mechanics，1962，14(1):136－157.

[2]SEGRE G，SILBERBERG A.Behaviour of macroscopic rigid spheres in poiseuille flow part 1.determination of local concentration by statistical analysis of particle passages through crossed light beams[J].Journal of Fluid Mechanics，1962，14(1):115 －135.

[3]JEFFREY R C，PEARSON J R A.Particle motion in laminar vertical tube flow[J].J.Fluid Mech，1965，22(4):721－735.

[4]KARNIS A，GOLDSMITH H L，MASON S G.The flow of suspensions through tubes:V.inertial effects[J].The Canadian Journal of Chemical Engineering，1966，44(4):181－193.

[5]MATAS J P，MORRIS J F，GUAZZELLI E.Inertial migration of rigid spherical particles in poiseuille flow[J].Journal of Fluid Mechanics，2004，515(1):171 －195.

[6]MATAS J P，MORRIS J F，GUAZZELLI E.Lateral forces on a sphere[J].Oil ＆ Gas Science and Technology，2004，59(1):59 －70.

[7]FENG J，HU H H，JOSEPH D D.Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid.part 2.couette and poiseuille flows[J].Journal of fluid mechanics，1994，277(271):271 －301.

[8]PAN T W，GLOWINSKI R.Direct simulation of the motion of neutrally buoyant circular cylinders in plane poiseuille flow[J].Journal of Computational Physics，2002，181(1):260－279.

[9]PAN T W，GLOWINSKI R.Direct simulation of the motion of neutrally buoyant balls in a three－dimensional poiseuille flow[J].Comptes Rendus Mécanique，2005，333(12):884－895.

[10]YU Z，PHAN-THIEN N，TANNER R I.Dynamic simulation of sphere motion in a vertical tube[J].Journal of Fluid Mechanics，2004(518):61 －93.

[11]SHAO X，YU Z，SUN B.Inertial migration of spherical particles in circular poiseuille flow at moderately high Reynolds numbers[J].Physics of Fluids，2008(20):103－307.

[12]CHUN B，LADD A J C.Inertial migration of neutrally buoyant particles in a square duct:an investigation of multiple equilibrium positions[J].Physics of Fluids，2006(18):317－320.

[13]YU Z，SHAO X.A direct-forcing fictitious domain method for particulate flows[J].Journal of computational physics，2007，227(1):292 －314.

[14]YU Z.A DLM/FD method for fluid/flexible-body interactions[J].Journal of Computational Physics，2005，207(1):1－27.