方管中中性懸浮顆粒慣性遷移現(xiàn)象的數(shù)值研究*

姚同林，余釗圣，邵雪明

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

0 引 言

顆粒懸浮流在自然界和工農(nóng)業(yè)生產(chǎn)中廣泛存在，如流化床、管道輸運(yùn)、細(xì)胞分離，以及更為常見(jiàn)的河道中泥沙沉降和輸運(yùn)等。盡管學(xué)術(shù)界對(duì)這種流動(dòng)已有了大量的實(shí)驗(yàn)、數(shù)值方面的研究，但由于顆粒與流體之間相互作用的復(fù)雜性，至今對(duì)其流動(dòng)特性以及作用機(jī)理仍缺乏足夠的認(rèn)識(shí)。

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

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

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

1 數(shù)值方法

1.1 物理模型

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

圖1 方管中一個(gè)顆粒的示意圖

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

對(duì)處于平衡態(tài)的方管中的流體，無(wú)量綱后的額外壓力梯度為:

1.2 虛擬區(qū)域法

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

基于虛擬區(qū)域法的假定，本研究在顆粒內(nèi)部引入虛擬體積力λ，假設(shè)顆粒的密度、體積、轉(zhuǎn)動(dòng)慣量、速度和角速度分別為 ρs，Vp，J，U 和 ωp，采用下述特征量來(lái)無(wú)量綱化控制方程:特征長(zhǎng)度Lc，特征速度Uc，特征時(shí)間Lc/Uc，特征密度 ρf，特征壓力，特征擬體力。

無(wú)量綱化之后的控制方程為:

式中:Ω—整個(gè)計(jì)算區(qū)域;P(t)—顆粒占據(jù)的區(qū)域;Ωf—流體區(qū)域，Ωf=ΩP(t);r—以顆粒中心為原點(diǎn)的位置矢徑。

無(wú)量綱化控制參數(shù)如下:

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

雷諾數(shù):Re= ρfUcLc/μ;

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

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

2 顆粒的慣性遷移分析

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

表1 單個(gè)顆粒算例的具體參數(shù)

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

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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

[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.