通道層流場中單顆粒慣性橫向遷移的數(shù)值研究

魏良俊 ，王昊利

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

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

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

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

1 顆粒兩相流基本方程

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

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

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

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

（2）動量方程：

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

1.2 離散相運(yùn)動方程

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

（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 力以及剪切梯度力等，其力學(xué)特征和相關(guān)數(shù)學(xué)表達(dá)分別如下：

（1） Magnus力

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

Magnus力表達(dá)式如下：

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

（2） Saffman力

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

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

（3） 剪切梯度力

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

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

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

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

3.1 通道模擬與邊界設(shè)置

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

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

矩形管當(dāng)量直徑公式為：

式中：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)格劃分與遷移計(jì)算

3.2.1 圓管

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

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

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

模擬計(jì)算前首先用ICEM 軟件對圓管進(jìn)行網(wǎng)格劃分，得到161 萬六面體結(jié)構(gòu)化網(wǎng)格，如圖3 所示。為研究顆粒橫向遷移情況，以圓管半徑R為特征長度對顆粒的位置信息進(jìn)行無量綱處理，假設(shè)顆粒在某時刻所處位置坐標(biāo)為 （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)流體及流向與圓管相同。模型的坐標(biāo)原點(diǎn)建立在方管的進(jìn)口截面中心處，坐標(biāo)系同圖1。顆粒釋放的初始位置位于進(jìn)口截面的z0=0、0.22L、0.33L處，如圖4 所示。圖4 中L為方管正方形截面邊長，z0為顆粒初始位置的z軸坐標(biāo)。

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

使用ICEM 軟件對方管進(jìn)行網(wǎng)格劃分，得到156萬個六面體結(jié)構(gòu)化網(wǎng)格，如圖5所示。為研究顆粒橫向遷移情況，以正方形截面邊長L為特征尺度對顆粒的位置信息進(jìn)行無量綱處理。根據(jù)文獻(xiàn)［13］報(bào)道的顆粒遷移數(shù)據(jù)以及本文大量仿真模擬結(jié)果，當(dāng)顆粒釋放位置在方管入口截面的中心對稱線上時，其在管道內(nèi)的遷移軌跡近似在同一平面內(nèi)。假設(shè)顆粒釋放位置為（0，0，z0），在某時刻顆粒所處的位置坐標(biāo)為（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 時，通過流體動力學(xué)軟件FLUENT（2021 R1）對圓管中的流體速度場進(jìn)行模擬，得到圓管內(nèi)流體流向最大截面的局部速度，如圖6 所示。從圖6 可以看到，速度場經(jīng)過了短暫的進(jìn)口段之后，進(jìn)入充分發(fā)展段。

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

4.1.2 圓管中顆粒的橫向遷移

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

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

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

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

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

圖10 給出了Re=100、不同直徑的顆粒在r0=0.8R處釋放時的遷移軌跡。從圖10 可以看到，當(dāng)d=0.8、1.0、1.2 mm 時，顆粒分別在r=0.660R、0.658R、0.657R處達(dá)到平衡。顯然，隨著粒子尺寸的增加，圓管內(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 時，通過流體動力學(xué)軟件FLUENT（2021 R1）對方管中的流體速度場進(jìn)行模擬，得到方管內(nèi)流體流向最大截面的局部速度，如圖11所示。從圖11可以看到，速度場經(jīng)過了短暫的進(jìn)口段之后，進(jìn)入充分發(fā)展段。

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

4.2.2 方管中顆粒的橫向遷移

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

圖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 時，顆粒在方管進(jìn)口截面中心即z0=0 處釋放的遷移軌跡。從圖12 可以看到，方管與圓管規(guī)律相似，均存在橫向位移現(xiàn)象，但數(shù)值可以忽略不計(jì)，即不具有明確的物理意義。

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

圖14 給出了不同雷諾數(shù)條件下粒徑d=1.0 mm 的顆粒在方管進(jìn)口截面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āng)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ù)量對計(jì)算結(jié)果的影響

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

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

圖17給出了網(wǎng)絡(luò)數(shù)量分別為125 萬、156 萬、187 萬時，顆粒在方管進(jìn)口截面z0=0.33L處釋放的遷移軌跡。從圖17 可以看到，3 種網(wǎng)格下顆粒遷移最終在z=0.253L～0.257L處達(dá)到平衡。顯然，不同網(wǎng)格方管計(jì)算結(jié)果偏差在合理范圍內(nèi)，說明方管計(jì)算結(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）對于圓管，當(dāng)雷諾數(shù)Re取100、125、150，顆粒直徑d采用0.8、1.0、1.2 mm，釋放位置r0=0.7R、0.8R時，顆粒慣性橫向遷移的平衡位置約在r=0.6R附近，且隨著粒子尺寸的增加，顆粒平衡位置將向管道軸心偏移；隨著Re的增加，顆粒平衡位置向管壁方向偏移。

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

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