微通道中顆粒慣性聚集特性的數(shù)值研究

石慧霞， 王企鯤

（上海理工大學(xué)能源與動力工程學(xué)院，上海 200093）

微通道中顆粒慣性聚集特性的數(shù)值研究

石慧霞， 王企鯤

（上海理工大學(xué)能源與動力工程學(xué)院，上海 200093）

根據(jù)運動相對性原理提出了一個描述單顆粒在無限長通道內(nèi)穩(wěn)定運動的準(zhǔn)定常計算模型，采用計算流體力學(xué)軟件研究了具有一定尺寸的球狀顆粒在方形截面微通道中所受到橫向升力的空間分布特征及其影響因素，并分析了其組成成分.研究結(jié)果表明，顆粒受到的橫向升力在通道y軸方向呈規(guī)則的空間分布，在數(shù)值上先正后負且存在唯一的零點，這即是顆粒的慣性聚集點；橫向升力主要由壓力橫向升力和剪切橫向升力組成，而壓力橫向升力是使顆粒產(chǎn)生慣性聚集現(xiàn)象的決定性因素.

慣性聚集；橫向升力；剛性球狀顆粒；低雷諾數(shù)流動；微通道流動；數(shù)值研究

當(dāng)細顆粒形成的均勻分布的流體以低雷諾數(shù)Re層流流入圓管時，經(jīng)過一段距離的流動后，這些顆粒會被穩(wěn)定地聚集在一個離管道中心距離為0.6倍半徑的同心圓環(huán)上，隨主流一起運動.這個運動顆粒圓環(huán)最初由Segre和Silberberg［1］通過實驗發(fā)現(xiàn)，常被稱為Segre-Silberberg圓環(huán)［2］.它表明顆粒在低Re的層流管內(nèi)運動時，除了受到沿主流方向的驅(qū)動力外，還受到垂直于主流方向的橫向升力的作用.驅(qū)動力會使顆粒加速運動，直到和周圍流體的速度相同，而橫向升力會引起顆粒產(chǎn)生垂直于主流方向的橫向遷移.這種橫向升力被認為是由于流場慣性力作用而產(chǎn)生的，故由其引發(fā)的聚集現(xiàn)象稱顆粒慣性聚集現(xiàn)象［2-3］.

目前，顆粒慣性聚集現(xiàn)象可以和微流控芯片技術(shù)相結(jié)合，研制出各種固液分離裝置，能被用于水的凈化和血細胞的篩選.這種固液分離裝置無需施加電、磁場或過濾膜等額外裝置，且結(jié)構(gòu)簡單，因而具有潛在的商業(yè)應(yīng)用價值［3-4］.Bhagat等［5］利用大高寬比的矩形截面微通道分離裝置，根據(jù)不同尺寸顆粒在同一Re下的聚集位置不同，成功地從尺寸為1.9～590μm的顆?；旌狭髦蟹蛛x出590μm的顆粒.高寬比為L／LC，L為通道高度，LC為通道寬度.

除直管外，彎管也常被用于分離裝置.Kuntaegowdanahalli等［6］利用顆粒慣性聚集原理，設(shè)計了一種螺旋形通道分離裝置，實現(xiàn)了尺寸為10，15，20μm這3種聚苯乙烯顆粒的連續(xù)分離，分離率達到90%. Carlo等［7］設(shè)計了一種非對稱正弦線彎管分離裝置，實現(xiàn)了不同尺寸聚二甲基硅氧烷（PDMS）顆粒、可變形硅油微粒和血液中血小板的分離，分離后的顆粒純度可達90%～100%，體積通量達到1 ml／min.

黃煒東等［8］利用慣性聚集原理，設(shè)計制作了一種具有不對稱彎管結(jié)構(gòu)的微流控芯片，并用它實現(xiàn)了對血液中血漿的分離.他們通過與現(xiàn)有的各種血漿分離方法相比較，發(fā)現(xiàn)采用顆粒慣性聚集原理制成的生物芯片能快速、有效地分離血漿，能使血液中血紅細胞的分離效率達到90%，同時基本不損害血細胞的活性.

盡管顆粒的慣性聚集原理在應(yīng)用層面上已作了大量研究，并擁有許多成功的案例，但是，其機理性研究，特別是關(guān)于對顆粒聚集起決定性作用的橫向升力的研究尚不完善，還有待深入［2-4，8］.

Asomolov［9］利用漸進匹配展開法給出了橫向升力的計算式，但該方法是將顆粒當(dāng)成點顆粒來處理的，忽略了顆粒自身尺寸對流場的影響，具有很大的局限性［2］.Matas等［10］通過研究Re=67～1 700情況下顆粒在Poiseuille流中的運動，揭示了隨著Re、顆粒尺寸的變化，顆粒慣性聚集位置的分布規(guī)律. Carlo等［11］通過數(shù)據(jù)擬合的方法得出了橫向升力在通道中心線附近與壁面附近的量級表達式，并發(fā)現(xiàn)旋轉(zhuǎn)對顆粒慣性聚集現(xiàn)象沒有顯著的影響.

本文建立了一個相對運動模型，利用計算流體力學(xué)（CFD）軟件［12-13］研究當(dāng)顆粒在方形截面通道y軸方向各個位置上的運動達到穩(wěn)定時，所受到的橫向升力沿通道y軸方向的空間分布規(guī)律，這為進一步深入研究橫向升力的力學(xué)特征及其對顆粒最終聚集位置的影響規(guī)律提供了理論指導(dǎo).

1 數(shù)值模擬

1.1 顆粒的相對運動模型

本文主要研究顆粒在方形通道中運動時，其在通道y軸各個位置處所受的橫向升力的變化規(guī)律.建立CFD數(shù)值計算模型時，如果以靜止的地面為參考系，該模型為非定常、動邊界問題.這種模型計算復(fù)雜、計算量極大，且不能針對性地揭示橫向升力沿通道y軸方向的空間分布.

鑒于本文研究的重點是顆粒在通道y軸方向的橫向升力，而不是準(zhǔn)確捕捉其在通道中的運動軌跡，因此，本文采用相對運動模型來進行準(zhǔn)定常的CFD計算.如圖1所示，坐標(biāo)系建立在顆粒上，并把顆粒置于需要獲得橫向升力的y軸位置處.當(dāng)顆粒運動到穩(wěn)定狀態(tài)時，其沿主流方向的速度為定值，此時它所受的垂直于主流方向的力即為橫向升力.采用相對運動模型能將非定常動邊界問題轉(zhuǎn)化為定常穩(wěn)態(tài)問題，這不僅減少了計算量，而且提高了計算精度.

圖1 顆粒的相對運動模型Fig.1 Relative motion model of particles

Carlo等［11］的研究表明，旋轉(zhuǎn)所產(chǎn)生的Magnus效應(yīng)對顆粒慣性聚集現(xiàn)象不會產(chǎn)生質(zhì)的改變，鑒于計算量的考慮，本文在計算時忽略了顆粒的旋轉(zhuǎn)，僅研究顆粒平動時的橫向升力分布規(guī)律.如圖1所示，在相對坐標(biāo)系中直徑為a的球狀顆粒即可被視為靜止不動，橫截面尺寸為H=50μm的方形通道四壁則以相對顆粒-Up的速度作勻速運動.在相對運動模型中，速度Up是個未知量，而在數(shù)值計算中Up要作為邊界條件來給定.因此，本文通過試湊的方法來得到速度Up，其具體過程如下：

a.給定顆粒的速度為U1進行數(shù)值計算，計算完成時，若顆粒在x軸方向所受到的驅(qū)動力為零，則認為顆粒在這一點的速度為Up=U1.如果驅(qū)動力小于零，則重新給定顆粒速度U2＜U1進行計算；反之，U2＞U1.在實際計算時，驅(qū)動力比升力低1～2個數(shù)量級就被認為驅(qū)動力為零.

b.步驟a計算完成后，若驅(qū)動力為零，則Up= U2；若驅(qū)動力小于零，則給定顆粒速度U3＜U2進行計算；反之，U3=（U1+U2）／2.

c.步驟b計算完成后，若驅(qū)動力為零，則Up= U3；若驅(qū)動力不為零，則重復(fù)步驟b中的方法進行計算，直到顆粒在x軸方向上的驅(qū)動力為零.此時，即可得到顆粒在通道y軸方向某一位置的速度Up.

1.2 數(shù)值模型的建立與其網(wǎng)格劃分

為保證顆粒所受的橫向升力是在流動為充分發(fā)展的管流流動基礎(chǔ)上數(shù)值計算得到的，本文將顆粒置于通道充分發(fā)展段區(qū)域.參照圓管層流進口段計算公式［14］近似得到方形通道層流進口段長度Le的計算公式為

式中，Le為層流起始段長度；H為通道水力直徑（即方形通道橫截面邊長）；，U為流體平均速度，ρ為流體密度，μ為流體動力粘性系數(shù).

在相對運動模型中，若方形通道的總長度大于Le，則可以忽略入口段對顆粒的影響.鑒于本文的研究內(nèi)容為顆粒在低Re（Re=20～80）下所受橫向升力的分布規(guī)律，通過計算選定通道總長為10 H.最終，計算模型的幾何尺寸設(shè)定為：通道長500μm，橫截面尺寸50μm×50μm.

考慮到計算精度，在方形通道的進口、出口區(qū)域生成結(jié)構(gòu)化網(wǎng)格，在顆粒周圍的區(qū)域生成非結(jié)構(gòu)化網(wǎng)格.由于顆粒附近區(qū)域的網(wǎng)格對橫向升力的影響較大，因此，增加該區(qū)域的網(wǎng)格數(shù)量.

為保證數(shù)值計算的準(zhǔn)確性，本文選擇靜止顆粒在通道中的繞流作為研究對象，以顆粒沿主流（在x軸方向）所受的驅(qū)動力作為衡量網(wǎng)格獨立性的驗證參數(shù).圖2為Re=0.5、直徑10μm的靜止顆粒在通道中心位置時所受驅(qū)動力系數(shù)CD隨網(wǎng)格數(shù)N的變化圖.驅(qū)動力計算式為

圖2 網(wǎng)格獨立性驗證結(jié)果Fig.2 Independence verification for grid

式中，P為應(yīng)力張量，計算式為P=-pI+P′，p為流體壓力，P′為偏應(yīng)力張量，I表示二階單位張量；i為主流方向（x軸方向）單位矢量；Σ為顆粒表面；n為法向單位矢量.

對應(yīng)的驅(qū)動力系數(shù)

式中，A為顆粒迎風(fēng)面積，A=πa2／4.

從圖2可以看出，當(dāng)網(wǎng)格數(shù)從60萬變到94萬時，驅(qū)動力系數(shù)CD的變化很小，表明網(wǎng)格數(shù)大于60萬時對計算結(jié)果影響很小，可認為60萬網(wǎng)格數(shù)是網(wǎng)格獨立解，因此，本文取60萬作為最終的計算網(wǎng)格數(shù).

1.3 數(shù)值解法及邊界條件

利用CFD軟件進行數(shù)值模擬，求解方法采用非耦合隱式求解.在離散化格式中，壓力速度耦合差值采用Simplec算法，壓力差值格式采用Standard算法，動量方程中對流項采用指數(shù)格式，粘性項采用中心差分格式.在邊界條件設(shè)置中，采用相對運動模型，顆粒靜止不動，通道壁面以速度Up沿來流反方向運動，通道進口速度為U-Up，出口采用自由出流，參考壓力設(shè)定在出口端.

2 結(jié)果及分析

2.1 橫向升力的變化曲線

為便于討論，定義了以下幾個無量綱參數(shù)：

a.橫向升力系數(shù)CFL.

式中，F(xiàn)L為數(shù)值計算得到的顆粒所受的橫向升力，計算式為

式中，y為顆粒距通道中心線的距離.

圖3為橫向升力系數(shù)CFL沿y+的分布圖.圖3（a）和圖3（b）表明，橫向升力系數(shù)CFL隨著y+的變化在通道y軸方向呈規(guī)則的空間分布.從通道中心為零開始，隨著y+的增大，橫向升力系數(shù)CFL正向增大，在y+=0.2附近達到最大值，之后開始下降，在y+=0.4～0.5附近達到零點，達到零點后隨著y+的增大，反向增大，且反向數(shù)值遠大于正向數(shù)值.

圖3 橫向升力系數(shù)分布圖Fig.3 Distribution of lift coefficient

本文的計算模型是將顆粒置于通道中心上半部分區(qū)域，如圖4所示，y為顆粒質(zhì)心到通道中心的距離.由圖4中顆粒在通道中的空間分布及圖3中橫向升力系數(shù)CFL沿通道y軸先正后負的規(guī)律性分布可知，CFL大于零，表明顆粒受到的橫向升力方向指向通道壁面；反之，指向通道中心.因此，顆粒在通道中運動時，在通道中心附近受到橫向升力的作用向通道壁面方向遷移，而在通道壁面附近在橫向升力的作用下向通道中心線方向遷移.

圖4 顆粒在通道中的分布Fig.4 Position of a particle in channel

橫向升力沿通道y軸方向分布，在顆粒從通道中心向壁面移動的過程中，其大小和方向都會發(fā)生改變.在這個過程中，便會出現(xiàn)零點位置.這個升力零點所在位置即是顆粒的慣性聚集點，顆粒在這個位置會發(fā)生慣性聚集.

圖5為顆粒尺寸a+=0.2時，系數(shù)CFL隨著Re的變化曲線圖.由圖5（a）和圖5（b）可知，在壁面附近區(qū)域（y+＞0.6）以及通道中心附近區(qū)域（y+＜0.4），橫向升力系數(shù)隨著Re的增大而逐漸減?。欢趹T性聚集點附近（0.4＜y+＜0.6），升力系數(shù)隨著Re的增大而增大.由于橫向升力在這個區(qū)域會出現(xiàn)零點，因此，在這個位置系數(shù)CFL變化很大.

2.2 聚集平衡位置

根據(jù)圖3（a）和圖3（b）中橫向升力系數(shù)CFL的變化曲線可知，曲線與橫坐標(biāo)的交點即為顆粒的慣性聚集點.因此，本文根據(jù)圖3中橫向升力沿通道y軸方向的變化曲線，通過插值得到橫向升力為零時顆粒在通道y軸方向的位置，該位置即可被認為是顆粒慣性聚集點，如圖6所示.

圖6（a）為顆粒慣性聚集位置隨顆粒尺寸a+的分布圖，該圖表明，隨著顆粒尺寸a+的增大，慣性聚集位置向通道中心位置移動，這與Matas等［10］及Carlo等［11］的研究結(jié)果一致.

圖6（b）為顆粒慣性聚集位置隨Re的變化圖.在Re=10～80的變化范圍內(nèi)，隨著Re的增大，顆粒在通道中的聚集位置變化趨勢不明顯.

2.3 流體壓力與粘性力對橫向升力的影響

顆粒在不可壓縮粘性流場中受到的力大致可以分為兩類：表面力和質(zhì)量力.在本文中，表面力主要表現(xiàn)為流體壓力和剪切應(yīng)力，質(zhì)量力主要表現(xiàn)為顆粒的重力，而本文不考慮顆粒的重力，因此，橫向升力只包括流體壓力和剪切應(yīng)力兩部分.

圖5 a＋＝0.2時，系數(shù)CFL隨Re的分布圖Fig.5 Coefficient CFLinfluenced by Re when a＋＝0.2

圖6 聚集平衡位置隨a＋，Re數(shù)的變化Fig.6 Focus position influenced by a＋and Re

將橫向升力中的流體壓力部分稱為壓力橫向升力，用符號FLP表示，計算式為

對應(yīng)的無量綱升力系數(shù)

將橫向升力中的剪切應(yīng)力部分稱為剪切橫向升力，用符號FLτ表示，計算式為

對應(yīng)的無量綱升力系數(shù)

圖7為Re=20時，壓力橫向升力系數(shù)CFP、剪切橫向升力系數(shù)CFτ隨y+的分布圖.圖7（a）和圖7（b）表明，在顆粒尺寸a+變化的過程中，壓力橫向升力系數(shù)CFP的分布規(guī)律與橫向力系數(shù)CFL相似.從通道中心為零開始，隨著y+的增大，壓力橫向升力系數(shù)CFP正向增大，在y+=0.2～0.3位置達到最大值，之后開始下降，在y+=0.5～0.6附近達到零點，達到零點后隨著y+的增大，反向增大，且反向數(shù)值遠大于正向數(shù)值.這表明壓力橫向升力在通道中心附近方向指向通道壁面，而在通道壁面附近方向指向通道中心.

當(dāng)顆粒尺寸較大時，剪切橫向升力系數(shù)CFτ隨著y+的改變，其數(shù)值幾乎都小于零，如圖7（b）所示.這表明隨y+的變化，剪切橫向升力的方向一直指向通道中心，而壓力橫向升力從通道中心到壁面先正后負且存在唯一的零點位置，這是產(chǎn)生顆粒慣性聚集現(xiàn)象的根本原因.

圖7 Re＝20時，系數(shù)CFP、系數(shù)CFτ隨y＋的分布Fig.7 Distribution of CFPand CFτwhen Re＝20

2.4 顆粒的速度變化曲線

顆粒在通道中隨流體運動時，受到流體驅(qū)動力的作用，以一定的速度沿主流方向運動.顆粒在通道中的速度主要取決于顆粒所在的y軸位置y+和Re.圖8為通過相對運動方法得到Re為20，80時顆粒在通道中的無量綱速度U+p分布圖.

圖8表明，顆粒在通道中的速度分布與流體的速度分布相似，在通道中心速度最大，從中心到通道壁面逐漸減小.在同一Re下，顆粒在通道同一位置的速度隨著a+的增大略有下降.從圖8（a）和圖8（b）中不難發(fā)現(xiàn)，在不同Re下，顆粒的速度變化趨勢相同，無量綱速度值也相近.

圖8 顆粒的速度分布圖Fig.8 Radial distribution of particles velocities

3 結(jié) 論

基于一個相對運動模型，通過CFD技術(shù)，數(shù)值研究了具有一定尺寸的顆粒在方形截面微通道中所受橫向升力的空間分布特征，并對橫向升力組成部分作了深入的研究，得到如下結(jié)論：

a.顆粒在通道中隨流體運動時，會受到具有規(guī)律性空間分布特征的橫向升力作用.顆粒在通道中心附近受到的橫向升力方向指向通道壁面，而在通道壁面附近方向指向通道中心.

b.橫向升力沿通道y軸分布并不均勻且存在唯一零點，這個升力零點即是顆粒的慣性聚集點.隨著顆粒尺寸a+的增大，慣性聚集點向通道中心移動.

c.橫向升力主要由壓力橫向升力和剪切橫向升力兩部分組成.壓力橫向升力的數(shù)值沿y+先正后負且存在零點，而剪切橫向升力均為負值.因此，顆粒慣性聚集現(xiàn)象完全取決于壓力橫向升力的作用.

［1］ Segre G，Silberberg A.Radial particle displacements in poiseuille flow of suspension［J］.Nature，1961，189（4760）：209-210.

［2］ Carlo D D.Inertial microfluidics［J］.Lab Chip，2009，9（21）：3038-3046.

［3］ 王企鯤，孫仁.管流中顆?！皯T性聚集”現(xiàn)象的研究進展及其在微流動中的應(yīng)用［J］.力學(xué)進展，42（6），2012：692-703.

［4］ 項南，朱曉璐，倪中華.慣性效應(yīng)在微流控芯片中的應(yīng)用［J］.化學(xué)進展，2011，23（9）：1945-1958.

［5］ Bhagat A A S，Kuntaegowdanahalli S S，Papautsky I. Inertial microfluidics for continuous particle filtration and extraction［J］.Microfluid Nanofluid，2009，7（2）：217-226.

［6］ Kuntaegowdanahalli SS，Bhagat A A S，Kumar G，et al.Inertial microfluidics for continuous particle separation in spiral microchannels［J］.Lab Chip，2009，9（20）：2973-2980.

［7］ Carlo D D，Edd J F，Irimia D，et al.Equilibrium separation and ltration of particles using differential inertial focusing［J］.Anal Chem，2008，80（6）：2204 -2211.

［8］ 黃煒東，張何，徐濤，等.基于慣性微流原理的微流控芯片用于血漿分離［J］.科學(xué)通報，2011，56（21）：1711 -1719.

［9］ Asmolov E S.The inertial lift on a spherical particle in a plane poiseuille flow at large channel reynolds number［J］.Fluid Mech，1999，381（4）：63-87.

［10］ Matas J，Morris J F，Guazzelli E.Inertial migration of rigid spherical particles in poiseuille flow［J］.Fluid Mech，2004，515（18）：171-195.

［11］ Carlo D D，Edd J F，Humphry K J，et al.Particle segregation and dynamics in confined flow［J］.Phys Rev Lett，2009，102（9）：094503（1-4）.

［12］ 姚征，陳康民.CFD通用軟件綜述［J］.上海理工大學(xué)學(xué)報，2002，24（2）：137-144.

［13］ 晁東海，郭雪巖.大顆粒氣固流化床內(nèi)兩相流動的CFD模擬［J］.上海理工大學(xué)學(xué)報，2010，32（4）：333 -339.

［14］ 歸柯庭，汪軍，王秋穎.工程流體力學(xué)［M］.北京：科學(xué)出版社，2003.

（編輯：石 瑛）

Numerical Investigation on Inertial Focus of Particles in Micro Channel

SHIHui-xia， WANGQi-kun
（School of Energy and Power Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China）

According to the principle of relative motion，a quasi-steady numerical model was proposed to describe the stable motion of a rigid spherical particle in an infinite straight channel. Based on the model and by using CFD（computational fluid dynamics）technique，the spatial features and influence mechanism of lateral lift were numerically investigated when a rigid spherical particle with certain size moving in a square cross-sectional micro channel，and the components of lift were also discussed.The results indicate that the lateral lift shows a regular spatial distribution along the y direction.Meanwhile the magnitude of the lift changes from positive to negative in the radial direction，with unique radial position of zero lift，which is just the equilibrium position for the inertial focus of particles.The lateral lift consists of two components，the pressure lift and shear lift，in which the pressure lift is the determinant factor to the inertial focus of particles.

inertial focus；lateral lift；rigid spherical particle；low Reynolds number flow；micro channel flow；numerical investigation

O 352

A

1007-6735（2013）04-0355-06

2012-12-28

國家教育部博士點青年基金資助項目（20113120120003）

石慧霞（1986-），女，碩士研究生.研究方向：微流體機械內(nèi)部流動.E-mail：shx1988＠foxmail.com

王企鯤（1978-），男，副教授.研究方向：低雷諾數(shù)粘性流體力學(xué).E-mail：wangqk＠usst.edu.cn