魏良俊 ，王昊利

1.鹽城工學院 機械工程學院，江蘇 鹽城 224051；2.金陵科技學院 機電工程學院，江蘇 南京 211169

21世紀60年代，Segré等［1］首先發(fā)現(xiàn)了均勻散布細顆粒的流體以層流形式進入直管，經(jīng)過長距離的流動后細顆粒會自發(fā)地進行徑向遷移，最終受慣性力的影響聚集在管道內(nèi)部的環(huán)形空間。這種現(xiàn)象稱為“慣性聚集”，其主要特點是無需施加任何外力，就可以在微通道內(nèi)對一定尺寸的顆粒進行聚集，從而實現(xiàn)對其分離與過濾。

為了解釋顆粒在管內(nèi)層流中遷移形成的慣性聚集現(xiàn)象，國內(nèi)外學者對顆粒在宏觀流場中運動所受作用力進行了大量的研究，提出顆粒在慣性層流場中受到Magnus 力［2］、Saffman 力［3］和剪切梯度力［4］等橫向力的共同作用。早期有學者運用“攝動法”等理論方法研究了彈性顆粒在低雷諾數(shù)的簡單剪切流場中的力學特性，發(fā)現(xiàn)彈性顆粒外形的變化也會導致顆粒產(chǎn)生橫向遷移［5-6］；Carlo［7-8］等研究表明，在不考慮顆粒自轉(zhuǎn)的情況下，“慣性聚集”現(xiàn)象依然存在，引發(fā)“慣性聚集”現(xiàn)象的橫向遷移是流場慣性力作用的結(jié)果；Zhao等［9］利用計算流體動力學（computational fluid dynamics，CFD）方法模擬了具有弱彈性的顆粒在二維通道內(nèi)的“慣性聚集”特點，表明顆粒的彈性變形會產(chǎn)生一個指向管道軸心的附加橫向升力，且越靠近壁面，附加升力就越大；Nakayama 等［10］通過實驗測量和數(shù)值模擬對圓管中懸浮顆粒的慣性聚集現(xiàn)象進行了研究，發(fā)現(xiàn)隨著雷諾數(shù)的增加存在3 個區(qū)間對應著不同的慣性聚集分布；Miura等［11］研究了雷諾數(shù)為100～1 200 的方形流道中懸浮球形顆粒的慣性遷移，發(fā)現(xiàn)雷諾數(shù)大于250 時方管內(nèi)部截面上存在8 個平衡位置，分別位于通道截面中心和角落，并且隨著雷諾數(shù)的增加，顆粒愈發(fā)向通道中心和拐角處聚集；Raffiee 等［12］使用三維數(shù)值模擬研究了黏彈性流體中粒子的升力分布，預測了不同參數(shù)下的平衡位置，分析了不同平衡點的穩(wěn)定性；Yu等［13］采用虛擬域方法對方形直流道內(nèi)中性粒子的聚焦平衡位置進行了研究，發(fā)現(xiàn)隨著流體彈性的增加，平衡位置連續(xù)出現(xiàn)在橫截面的中線、對角線、拐角和流道的中心；全列［14］采用了格子玻爾茲曼方法，對中性懸浮球形顆粒在不同形狀管道中的輸運過程進行數(shù)值模擬，發(fā)現(xiàn)隨著雷諾數(shù)的增加多顆粒成鏈現(xiàn)象更加明顯，在軸對稱的圓管和矩形管中顆粒慣性遷移軌跡在橫截面上的投影分別是一條直線和一條曲線；高增鋒［15］利用FLUENT 對二維平面通道中單顆粒的慣性遷移進行了模擬，指出了黏性流體中顆粒慣性遷移現(xiàn)象的根本原因和顆粒所受慣性升力沿顆粒的橫向位置呈現(xiàn)規(guī)律性的空間分布。

為研究顆粒慣性聚集機制，本文采用計算流體動力學方法對圓管和方管內(nèi)隨流體運動的單個顆粒的橫向遷移進行模擬研究。將顆粒放在進口截面的若干不同位置并釋放，采用離散顆粒模型（discrete particle model，DPM）計算顆粒在流場中的遷移軌跡及橫向遷移后所達到的平衡位置。

1 顆粒兩相流基本方程

對于顆粒兩相流，分別考慮連續(xù)相（流體）控制方程與離散相（顆粒）運動方程。

1.1 連續(xù)相控制方程

連續(xù)相控制方程考慮不可壓縮的牛頓流體，滿足連續(xù)性方程和動量方程，分別如下：

（1）連續(xù)性方程：

（2）動量方程：

式中：?為Hamilton 算子；uf為流體流動速度，m/s；t為流體在管道中的流動時間，s；p為流體的壓強，Pa；ρf為流體密度，kg/m3；v為流體運動黏度，m2/s。

1.2 離散相運動方程

對于離散相模型，基于拉格朗日坐標系，顆粒從進口點開始跟蹤，直到離開流場或者滿足某種積分極限準則為止。離散相方程包括：

（1）顆粒軌跡方程：

（2）顆粒動量方程：

式中：Xp為顆粒的質(zhì)心位置，m；up為顆粒的速度，m/s；mp為顆粒的質(zhì)量，kg；g為重力加速度，m/s2；Fd為顆粒受到的繞流阻力，N；FVM為顆粒受到的附加質(zhì)量力，N；FB為Basset 力，N；FP為壓力梯度力，N；FML為Magnus 力，N；FS為Saffman 力，N；FX為其余外力的綜合，N。

2 顆粒橫向受力分析

根據(jù)已有的研究，顆粒在層流場中的橫向受力主要包括Magnus 力、Saffman 力以及剪切梯度力等，其力學特征和相關(guān)數(shù)學表達分別如下：

（1） Magnus力

當一個旋轉(zhuǎn)物體的角速度矢量與位移速度矢量不重合時，在與角速度矢量和位移速度矢量組成的平面相垂直的方向上產(chǎn)生一個橫向力，這個橫向力就是Magnus 力。在Magnus 力的作用下，物體位移軌跡會發(fā)生偏轉(zhuǎn)，這種現(xiàn)象稱作馬格努斯效應。

Magnus力表達式如下：

式中：d為顆粒直徑，m；ω為流體旋轉(zhuǎn)的角速度，rad/s。

（2） Saffman力

流過顆粒表面的流體具有不同的速度，在顆粒表面會產(chǎn)生一個由低速流體側(cè)指向高速流體側(cè)的橫向作用力，這個橫向作用力就是Saffman力，表達式如下：

式中：K為一常數(shù)，通常取6.46；μ為流體動力學黏度，Pa·s。

（3） 剪切梯度力

根據(jù)伯努利原理，當流體對顆粒具有相對運動速度時，在垂直于顆粒運動方向上會產(chǎn)生一個力，這個力就是剪切梯度力，表達式如下［16］：

式中：G為剪切梯度力，N；ΔP為顆粒上下表面的壓力差，Pa；ΔH為顆粒上下表面的距離，m。

3 顆粒橫向遷移數(shù)值模擬

為研究顆粒橫向遷移情況，本文以中性懸浮單顆粒在液相層流通道內(nèi)的穩(wěn)定運動為研究對象，忽略質(zhì)量力和顆粒對流場的影響，采用“顆粒點源”直接數(shù)值模擬方法進行通道層流場中顆粒的慣性遷移研究。該方法將顆粒當作一個質(zhì)點，對液相采用歐拉法建立控制方程，用拉格朗日方法模擬顆粒在流場中的運動；在流場計算基礎上，對顆粒引入Saffman力和剪切梯度力進行計算和軌跡追蹤，暫未考慮顆粒旋轉(zhuǎn)所產(chǎn)生的Magnus力的影響，主要原因是顆粒旋轉(zhuǎn)作用效果有限。

3.1 通道模擬與邊界設置

基于現(xiàn)有條件，采用計算流體動力學軟件FLUENT（2021 R1）對圓管和方管在層流場中的顆粒慣性遷移開展數(shù)值模擬研究。邊界條件為：進口速度條件設置為截面平均速度，具體數(shù)值根據(jù)設定的Re，按照式（8）反算計算求得；出口為壓力出口；壁面為無滑移邊界條件。在計算過程中，壓力、速度的耦合采用SIMPLE算法計算。

式中：U為流經(jīng)通道截面流體的平均速度，m/s；Re為雷諾數(shù)；l為管道特征長度，對于圓管為半徑，方管取當量直徑，m。

矩形管當量直徑公式為：

式中：a、b分別為矩形管矩形截面中兩條邊的長度，m。

表1 給出了不同雷諾數(shù)下，流經(jīng)圓管和方管截面流體的平均速度。

表1 不同雷諾數(shù)下通道截面流體的平均速度Table 1 Average velocity of channel section fluid at different Re m/s

3.2 網(wǎng)格劃分與遷移計算

3.2.1 圓管

圓管通道模型中的圓管直徑為9.6 mm，長度為1.0 m，管內(nèi)流體為水。模型的坐標原點建在圓管進口截面中心處，如圖1 所示，其中X軸為流體的流動方向，Y和Z為橫向坐標軸。顆粒釋放的初始位置位于進口截面的r0=0、0.7R、0.8R的點，如圖2 所示。圖2 中R為圓管截面的半徑，r0為顆粒初始位置的半徑。

圖1 部分圓管模型及顆粒橫向受力示意圖Fig. 1 Schematic diagram of partial circular tube model and lateral forces on particles

圖2 顆粒在圓管進口處釋放時的位置示意圖Fig. 2 Schematic diagram of particle release position at the entrance of circular tube

模擬計算前首先用ICEM 軟件對圓管進行網(wǎng)格劃分，得到161 萬六面體結(jié)構(gòu)化網(wǎng)格，如圖3 所示。為研究顆粒橫向遷移情況，以圓管半徑R為特征長度對顆粒的位置信息進行無量綱處理，假設顆粒在某時刻所處位置坐標為 （x1，y1，z1），則顆粒的無量綱流向與徑向的位移x′、r′分別為：

圖3 圓管局部網(wǎng)格圖Fig. 3 Local grid diagram of circular tube

3.2.2 方管

方管通道模型中的方管邊長為9.1 mm，長為1.0 m，管內(nèi)流體及流向與圓管相同。模型的坐標原點建立在方管的進口截面中心處，坐標系同圖1。顆粒釋放的初始位置位于進口截面的z0=0、0.22L、0.33L處，如圖4 所示。圖4 中L為方管正方形截面邊長，z0為顆粒初始位置的z軸坐標。

圖4 顆粒在方管進口處釋放時的位置示意圖Fig. 4 Schematic diagram of particle release position at the entrance of square pipe

使用ICEM 軟件對方管進行網(wǎng)格劃分，得到156萬個六面體結(jié)構(gòu)化網(wǎng)格，如圖5所示。為研究顆粒橫向遷移情況，以正方形截面邊長L為特征尺度對顆粒的位置信息進行無量綱處理。根據(jù)文獻［13］報道的顆粒遷移數(shù)據(jù)以及本文大量仿真模擬結(jié)果，當顆粒釋放位置在方管入口截面的中心對稱線上時，其在管道內(nèi)的遷移軌跡近似在同一平面內(nèi)。假設顆粒釋放位置為（0，0，z0），在某時刻顆粒所處的位置坐標為（x2，y2，z2），則顆粒的無量綱流向與徑向位移x″、z″分別為：

圖5 方管局部網(wǎng)格圖Fig. 5 Local grid diagram of square tube

4 模擬與分析

4.1 圓管

4.1.1 圓管內(nèi)速度場的模擬

為研究圓管內(nèi)顆粒橫向遷移規(guī)律，對直徑為9.6 mm、長度為1.0 m、管內(nèi)流體為水、具有161萬個六面體結(jié)構(gòu)化網(wǎng)格的圓管，在未加入顆粒、雷諾數(shù)Re=100 時，通過流體動力學軟件FLUENT（2021 R1）對圓管中的流體速度場進行模擬，得到圓管內(nèi)流體流向最大截面的局部速度，如圖6 所示。從圖6 可以看到，速度場經(jīng)過了短暫的進口段之后，進入充分發(fā)展段。

圖6 圓管內(nèi)流向最大截面的局部速度云圖Fig. 6 Local velocity nephogram of the maximum cross section of circular tube flow direction

4.1.2 圓管中顆粒的橫向遷移

在圖6 速度場基礎上，采用離散顆粒相計算方法，對粒徑d=0.8、1.0、1.2 mm 的中性球形顆粒，在雷諾數(shù)Re=100、125、150 條件下，在進口截面r0=0、0.7R、0.8R的3個初始位置釋放的遷移過程分別進行計算，得到圓管內(nèi)流場中的顆粒遷移軌跡及橫向遷移后所達到的平衡位置，如圖7～圖10 所示。計算時僅考慮流體對顆粒的單向作用，不考慮顆粒對流場的反作用。

圖7 顆粒在r0=0處釋放的遷移軌跡Fig. 7 Migration trajectory of particles release when r0=0

圖8 顆粒在r0=0.7R、0.8R處釋放的遷移軌跡Fig. 8 Migration trajectory of particles release when r0=0.7R and 0.8R

圖9 不同Re時顆粒在圓管內(nèi)的遷移軌跡Fig. 9 Migration trajectory of particles with different Re in a circular tube

圖10 不同粒徑時顆粒在圓管內(nèi)的遷移軌跡Fig. 10 Migration trajectory of particles with different sizes in a circular tube

圖7 給出了Re=100、d=0.8 mm，Re=150、d=1.0 mm，Re=100、d=1.0 mm 時，顆粒在圓管進口截面中心即r0=0 處釋放的遷移軌跡。從圖7 可以看到，雖然存在橫向遷移現(xiàn)象，但3種情況下顆粒最大橫向遷移距離僅為0.008R，可以忽略不計。產(chǎn)生偏移的主要原因可能是由于數(shù)值波動，不具有明確的物理意義。

圖8 給出了直徑d=1.0 mm 的顆粒在圓管進口截面r0=0.7R、0.8R處釋放時的遷移軌跡。從圖8 可以看到，顆粒遷移過程經(jīng)歷了初始和平衡兩個主要階段，在初始階段主要表現(xiàn)為顆粒跟隨主流向下游運動的同時發(fā)生橫向遷移；經(jīng)過短暫的運動后，顆粒的橫向遷移達到平衡階段，其中，從r0=0.7R處釋放的顆粒在r=0.57R的半徑處達到平衡，從r0=0.8R處釋放的顆粒在r=0.65R的半徑處達到平衡。

圖8 表明，顆粒在圓管內(nèi)被層流向下游驅(qū)動時，會自發(fā)地進行徑向遷移，最終會穩(wěn)定地停留在離管道軸心約0.6倍半徑的位置上并向下游繼續(xù)運動。該模擬結(jié)果與Segre 等［1］的實驗觀測結(jié)果基本吻合。

圖9 給出了粒徑d=1.0 mm 的顆粒在r0=0.8R處釋放時在不同雷諾數(shù)下的遷移軌跡。從圖9可以看到，當Re=100、125、150 時，顆粒分別在r=0.658R、0.665R、0.672R處達到平衡。顯然，隨著Re的增加，圓管內(nèi)顆粒慣性遷移的平衡位置向管壁方向偏移，而且從r0=0.7R處釋放的顆粒也具有相似的規(guī)律。

圖10 給出了Re=100、不同直徑的顆粒在r0=0.8R處釋放時的遷移軌跡。從圖10 可以看到，當d=0.8、1.0、1.2 mm 時，顆粒分別在r=0.660R、0.658R、0.657R處達到平衡。顯然，隨著粒子尺寸的增加，圓管內(nèi)顆粒慣性遷移的平衡位置將向管道軸心偏移，而且從r0=0.7R處釋放的顆粒也具有相似規(guī)律。

4.2 方管

4.2.1 方管內(nèi)速度場的模擬

為研究方管內(nèi)顆粒橫向遷移規(guī)律，對邊長為9.1 mm、長度為1.0 m、管內(nèi)流體為水、具有156萬個六面體結(jié)構(gòu)化網(wǎng)格的方管，在未加入顆粒、雷諾數(shù)Re=100 時，通過流體動力學軟件FLUENT（2021 R1）對方管中的流體速度場進行模擬，得到方管內(nèi)流體流向最大截面的局部速度，如圖11所示。從圖11可以看到，速度場經(jīng)過了短暫的進口段之后，進入充分發(fā)展段。

圖11 部分方管中心截面的局部速度云圖FIg. 11 Local velocity nephogram of central section of partial square tube

4.2.2 方管中顆粒的橫向遷移

在圖11 速度場基礎上，采用離散顆粒相計算方法，對粒徑d=0.8、1.0、1.2 mm 的中性球形顆粒，在雷諾數(shù)Re=100、125、150 條件下，在進口截面z0=0、0.22L、0.33L的3 個初始位置釋放的遷移過程分別進行計算，得到方管內(nèi)流場中的顆粒遷移軌跡及橫向遷移后所達到的平衡位置，如圖12～圖15所示。計算時僅考慮流體對顆粒的單向作用，不考慮顆粒對流場的反作用。

圖12 顆粒在z0=0處釋放的遷移軌跡Fig. 12 Migration trajectory of particles released at z0=0

圖13 顆粒在z0=0.22L、0.33L處釋放的遷移軌跡Fig. 13 Migration trajectory of particles released at z0=0.22L、0.33L

圖14 不同Re時顆粒在方管內(nèi)的遷移軌跡Fig. 14 Migration trajectory of particles with different Re in a square tube

圖15 不同粒徑時顆粒在方管內(nèi)的遷移軌跡Fig. 15 Migration trajectory of particles with different sizes in a square tube

圖12 給出了Re=100、d=1.0 mm，Re=150、d=1.0 mm，Re=100、d=0.8 mm 時，顆粒在方管進口截面中心即z0=0 處釋放的遷移軌跡。從圖12 可以看到，方管與圓管規(guī)律相似，均存在橫向位移現(xiàn)象，但數(shù)值可以忽略不計，即不具有明確的物理意義。

圖13 給出了顆粒在方管進口截面z0=0.22L、0.33L處釋放時的遷移軌跡。從圖13 可以看到，當釋放位置為z0=0.22L時，顆粒遷移的平衡位置約在z=0.17L處；當釋放位置為z0=0.33L時，顆粒遷移的平衡位置約在z=0.26L處。顯然，釋放位置越靠近管道軸心時，顆粒的橫向遷移距離越短，也越容易達到平衡位置。

圖14 給出了不同雷諾數(shù)條件下粒徑d=1.0 mm 的顆粒在方管進口截面z0=0.33L處釋放時的遷移軌跡。從圖14可以看到，Re=100、125、150時顆粒遷移平衡位置分別位于z=0.256L、0.255L、0.254L處。顯然，對于方管而言，在釋放位置等其他條件相同的情況下，Re不同，顆粒最終的平衡位置亦有所不同；隨著Re的增加，顆粒的平衡位置會更加偏向通道中心。顆粒從z0=0.22L處釋放時也具有相似的規(guī)律。

圖15給出了不同粒徑條件下，Re=100、顆粒在z0=0.33L處釋放時的遷移軌跡。從圖15 可以看到，當d=0.8、1.0、1.2 mm 時，顆粒遷移的平衡位置分別位于z=0.260L、0.256L、0.248L處。顯然，對于方管而言，在釋放位置等其他條件相同時，粒徑大小不同，顆粒遷移的平衡位置也不同，表現(xiàn)為小顆粒的平衡位置更加靠近壁面。顆粒從z0=0.22L處釋放時也具有相似的規(guī)律。

5 網(wǎng)格數(shù)量對計算結(jié)果的影響

圖16 給出了網(wǎng)絡數(shù)量分別為141 萬、161 萬、211 萬時，顆粒在圓管進口截面r0=0.8R處釋放的遷移軌跡。從圖16 可以看到，初始階段3 個網(wǎng)格的遷移軌跡完全吻合；經(jīng)過短暫遷移后，顆粒由初始r0=0.8R分別遷移至r=0.65R、0.64R、0.63R處達到平衡，相鄰兩種網(wǎng)格的計算偏差約為0.01R。顯然，計算結(jié)果未隨著網(wǎng)格數(shù)量增加產(chǎn)生明顯變化，即計算結(jié)果與網(wǎng)格數(shù)量無明顯關(guān)聯(lián)，計算所得的顆粒橫向遷移規(guī)律具有可信度。

圖16 不同網(wǎng)格數(shù)的圓管對顆粒遷移結(jié)果的影響Fig. 16 Influence on particles release in different grid numbers of circular tube

圖17給出了網(wǎng)絡數(shù)量分別為125 萬、156 萬、187 萬時，顆粒在方管進口截面z0=0.33L處釋放的遷移軌跡。從圖17 可以看到，3 種網(wǎng)格下顆粒遷移最終在z=0.253L～0.257L處達到平衡。顯然，不同網(wǎng)格方管計算結(jié)果偏差在合理范圍內(nèi)，說明方管計算結(jié)果與網(wǎng)格數(shù)量無顯著相關(guān)性。

圖17 不同網(wǎng)格數(shù)的方管對顆粒遷移結(jié)果的影響Fig. 17 Influence on particles release in different grid numbers of square tube

6 結(jié)論

通過對通道層流場中單顆粒橫向慣性遷移的數(shù)值研究，得出如下結(jié)論：

（1）對于圓管，當雷諾數(shù)Re取100、125、150，顆粒直徑d采用0.8、1.0、1.2 mm，釋放位置r0=0.7R、0.8R時，顆粒慣性橫向遷移的平衡位置約在r=0.6R附近，且隨著粒子尺寸的增加，顆粒平衡位置將向管道軸心偏移；隨著Re的增加，顆粒平衡位置向管壁方向偏移。

（2）對于方管，當雷諾數(shù)Re取100、125、150，顆粒直徑d采用0.8、1.0、1.2 mm，釋放位置z0=0.22L時，顆粒橫向遷移的平衡位置約在z=0.170L附近；當釋放位置z0=0.33L時，顆粒橫向遷移的平衡位置約在0.255L附近。顯然，釋放位置越靠近通道軸心，顆粒的橫向遷移距離越短；隨著Re的增加，顆粒的平衡位置會更加偏向管道中心。此外，小顆粒的平衡位置更靠近壁面。

（3）網(wǎng)格數(shù)量對計算結(jié)果偏差的影響均在合理范圍內(nèi)，計算結(jié)果具有可信度。