玻璃質(zhì)微通道流動(dòng)阻力特性的數(shù)值模擬

周萍，陳卓，徐則林，莫景文

(中南大學(xué) 能源科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙，410083)

玻璃因其特殊的物理化學(xué)性能，廣泛應(yīng)用于封裝、射頻微機(jī)電以及生物芯片基座等微流控系統(tǒng)中，因此，基于玻璃材質(zhì)的微加工工藝也就成為開(kāi)發(fā)微流控系統(tǒng)的一個(gè)重要技術(shù)問(wèn)題。玻璃微通道的刻蝕方法主要有機(jī)械刻蝕、干法刻蝕(等離子刻蝕)與濕法刻蝕三大類(lèi)。機(jī)械刻蝕包括傳統(tǒng)鉆刻、超聲鉆刻、電化學(xué)加工和噴粉法等，其缺點(diǎn)是無(wú)法得到較光滑的刻蝕表面。干法刻蝕的刻蝕成本高，刻蝕速率低，刻蝕面的表面粗糙度也較高[1]，因此，目前使用最廣泛的仍是濕法刻蝕。濕法刻蝕是利用氫氟酸對(duì)玻璃的腐蝕性對(duì)控制面的玻璃進(jìn)行腐蝕。由于玻璃各向同性的特征，刻蝕一般無(wú)法得到嚴(yán)格的矩形微通道截面形狀，而是側(cè)壁為圓弧形的通道[2-3]，本文將這種截面形狀的微通道稱(chēng)為圓角形微通道。同時(shí)，由于刻蝕的橫向速度遠(yuǎn)遠(yuǎn)大于縱向速度，因此，用一般濕法刻蝕方法得到的微通道截面大多寬度大于深度。雖然許多研究用復(fù)雜工藝可以得到玻璃濕法刻蝕的較理想矩形通道截面[4-5]，但由于經(jīng)濟(jì)性較差，目前尚無(wú)法用于批量生產(chǎn)。另一方面，在微流控系統(tǒng)設(shè)計(jì)中，阻力是非常重要的1個(gè)參數(shù)，其大小在很大程度上取決于微通道的截面形狀及其內(nèi)部流動(dòng)特征(層流、湍流以及過(guò)渡流)，常用范寧阻力系數(shù)來(lái)表示，即

式中：Δp為微通道壓降；Dh為截面特征尺度；L為微通道長(zhǎng)度；ρ為流體密度；u為流體流速。

范寧阻力系數(shù)(f)與雷諾數(shù)(Re)的乘積稱(chēng)為Poiseuille數(shù)(簡(jiǎn)稱(chēng) Po數(shù)，Po=f·Re)。對(duì)于圓管層流，Po是1個(gè)常數(shù)，其值為16；對(duì)于矩形通道，Po是1個(gè)與微通道截面形狀因子有關(guān)的量，Shah等[6]提出了在低雷諾數(shù)下Po的經(jīng)驗(yàn)公式，即：

微通道中流體的壓降受流動(dòng)狀態(tài)、流體性質(zhì)、通道截面形狀等眾多參數(shù)的影響，過(guò)去對(duì)該問(wèn)題的研究主要集中在幾個(gè)方面：

(1) 流體性質(zhì)對(duì)微通道流動(dòng)特性的影響。Wu等[7-8]對(duì)不可壓縮和可壓縮氮?dú)?、氫氣和氬氣的流?dòng)進(jìn)行了研究；Yang等[9]對(duì)不可壓縮空氣的流動(dòng)進(jìn)行了研究；Judy等[10]研究了異丙醇、水和甲醇的流動(dòng)；Yang等[9,11]研究了R134a液體和R114液的流動(dòng)。

(2) 截面形狀對(duì)微通道流動(dòng)特性的影響。Celata等[11]對(duì)圓形截面、馬璨等[12]對(duì)矩形截面進(jìn)行了研究，蔣炳炎等[13]對(duì)三角形和梯形截面微通道內(nèi)流動(dòng)特征進(jìn)行了研究。

(3) 管壁對(duì)微通道流動(dòng)特性的影響。王瑋等[14-16]通過(guò)研究發(fā)現(xiàn)，粗糙管壁會(huì)增加層流區(qū)的阻力系數(shù)以及引起層流向湍流提前轉(zhuǎn)變。

(4) 流動(dòng)狀態(tài)對(duì)微通道流動(dòng)特性的影響。

不同研究者分別對(duì)層流、過(guò)渡區(qū)流動(dòng)和湍流流動(dòng)[17]以及微通道內(nèi)的兩相流[18]，滑移流動(dòng)[19]等進(jìn)行了研究。已有研究工作主要是針對(duì)圓形、矩形、三角形、梯形等規(guī)則截面形狀的微通道而進(jìn)行的。對(duì)于實(shí)際玻璃濕法刻蝕得到截面形狀的微通道阻力特性的研究尚未見(jiàn)相關(guān)文獻(xiàn)報(bào)道。微通道流動(dòng)特性的研究方法主要有實(shí)驗(yàn)研究與數(shù)值模擬兩大類(lèi)。實(shí)驗(yàn)研究方法由于存在成本高、周期長(zhǎng)以及參數(shù)檢測(cè)困難等問(wèn)題而使其應(yīng)用受到很大的限制，數(shù)值模擬方法則由于具有成本低、周期短以及模擬能力強(qiáng)的特點(diǎn)，近些年在微通道流動(dòng)的研究領(lǐng)域得到了廣泛應(yīng)用[20-21]。為此，本文作者采用數(shù)值模擬的方法，以水為介質(zhì)，研究低雷諾數(shù)下圓角型微通道的截面形狀對(duì)阻力系數(shù)的影響規(guī)律，并擬合出Po準(zhǔn)數(shù)與界面形狀因子α之間的關(guān)聯(lián)式，為微通道流動(dòng)特性的研究及其工程應(yīng)用提供理論依據(jù)。

1 玻璃質(zhì)微通道幾何結(jié)構(gòu)及其網(wǎng)格劃分

根據(jù)濕法刻蝕的特點(diǎn)，數(shù)值模擬中對(duì)微通道截面形狀進(jìn)行適當(dāng)簡(jiǎn)化處理。設(shè)兩側(cè)壁為1/4圓弧面，底面為平面(見(jiàn)圖1)，且截面最大寬度D大于高度H。定義高H(亦即圓弧所在圓的半徑)與底部直段長(zhǎng)W的比值為截面形狀因子α，

結(jié)合作者前期工作，取微通道的寬度 D為 200 μm，長(zhǎng)度L為3 000 μm?？紤]W與H為8種不同尺寸即 8種不同的形狀因子，其具體的幾何參數(shù)見(jiàn)表1。

圖1 圓角形微通道幾何形狀示意圖Fig.1 Schematic diagram of geometry for semi-circle microchannel

表1 不同形狀因子α的圓角形微通道幾何參數(shù)Table 1 Geometric parameters of microchannels with different aspect ratios μm

對(duì)微通道的解析區(qū)域采用結(jié)構(gòu)化網(wǎng)格劃分，且對(duì)壁面處的網(wǎng)格加密以更好地處理邊界層流動(dòng)現(xiàn)象，使求解更精確。根據(jù)微通道截面幾何尺寸的不同，所劃分的總網(wǎng)格數(shù)為5.0×105～7.0×105，并對(duì)網(wǎng)格獨(dú)立性進(jìn)行驗(yàn)證。

2 數(shù)學(xué)模型及其邊界條件

2.1 數(shù)學(xué)模型

在對(duì)微通道流動(dòng)過(guò)程進(jìn)行數(shù)值模擬時(shí)，作出以下假設(shè)：(1) 流體為牛頓流體；(2) 流動(dòng)過(guò)程為等溫穩(wěn)態(tài)過(guò)程；(3) 流體為不可壓縮流體，且為層流流動(dòng)(0.01＜Re＜40)；(4) 流體在壁面附近為無(wú)滑移流動(dòng)；(5) 對(duì)于分子自由程小于10-3μm的流體，當(dāng)微通道特征尺度大于1 μm時(shí)，均能滿(mǎn)足克努曾準(zhǔn)數(shù)Kn＜10-3，流動(dòng)不受稀薄效應(yīng)影響，可以使用連續(xù)介質(zhì)假設(shè)[22]。所研究的微通道特征尺度均在30 μm以上，故其內(nèi)部流體流動(dòng)可視為連續(xù)流動(dòng)。

基于上述假設(shè)，描述微通道內(nèi)流體流動(dòng)過(guò)程的控制方程如下。

連續(xù)性方程：

動(dòng)量方程：

其中：ui為i方向上的速度分量，且i和j 均為1，2和3，代表x，y和 z 3個(gè)坐標(biāo)方向；p為靜壓；ρ和μ分別為流體的密度和動(dòng)力黏度。

2.2 邊界條件

微通道中的流動(dòng)介質(zhì)為水，溫度為定常值(293 K)，其密度 ρ=998.2 kg/m3，動(dòng)力黏度 μ=1.004×10-3Pa·s；流體入口為速度入口邊界條件，壁面為無(wú)滑移邊界條件；出口為壓力出口。

3 結(jié)果與討論

基于 FLUENT6.3軟件平臺(tái)，以圓形截面微通道內(nèi)流體阻力的理論計(jì)算公式為基準(zhǔn)，對(duì)微通道內(nèi)流體阻力數(shù)值模擬的模型與方法的有效性進(jìn)行檢驗(yàn)。在此基礎(chǔ)上，對(duì)截面形狀因子對(duì)圓角形微通道內(nèi)流動(dòng)阻力的影響規(guī)律進(jìn)行數(shù)值模擬研究。

3.1 雷諾數(shù)和截面形狀因子對(duì)微通道壓降的影響

根據(jù)文獻(xiàn)[23]和[24]中的研究經(jīng)驗(yàn)，對(duì) 0.01≤Re≤40時(shí)水在圓角形微通道內(nèi)的流動(dòng)進(jìn)行研究。針對(duì)微通道截面最大寬度D大于高度H的特點(diǎn)，取微通道的高度H為雷諾數(shù)的特征尺度[22]，且雷諾數(shù)的定義式為：

式中：um為微通道截面平均流速?；诓煌字Z數(shù)與截面形狀因子的微通道流動(dòng)過(guò)程的數(shù)值模擬結(jié)果，為排除入口和出口效應(yīng)的影響，由管流充分發(fā)展段計(jì)算出單位長(zhǎng)度壓降，并分別繪出微通道單位長(zhǎng)度壓降隨雷諾數(shù)、截面形狀因子變化的曲線，見(jiàn)圖2和圖3。

由圖2可知：不同截面形狀因子下圓角形微通道單位長(zhǎng)度壓降均隨雷諾數(shù)的增加而增加，且呈線性變化，這與常規(guī)層流流動(dòng)阻力的規(guī)律相同。圖3表明：截面形狀因子對(duì)玻璃質(zhì)微通道流動(dòng)阻力有重要影響；在微通道截面的最大寬度D一定的情況下，不同雷諾數(shù)下微通道單位長(zhǎng)度壓降均隨著截面形狀因子的增加而增加，且呈現(xiàn)出非線性變化的趨勢(shì)。

圖2 定寬度時(shí)雷諾數(shù)Re對(duì)微通道單位長(zhǎng)度壓降p的影響Fig.2 Effect of Re on micro-channel pressure drop at a constant width

圖3 定寬度時(shí)微通道形狀因子α對(duì)單位長(zhǎng)度壓降p的影響Fig.3 Effect of micro-channel shape factor on pressure drop at a constant width

3.2 雷諾數(shù)和截面形狀因子對(duì)微通道Po數(shù)的影響

由式(1)與(6)可得：

研究結(jié)果表明：當(dāng)截面形狀因子一定時(shí)，不同雷諾數(shù)下的Po重合在一起(如圖4所示，7個(gè)雷諾數(shù)條件下的數(shù)據(jù)點(diǎn)重合)，說(shuō)明Po并不隨Re的變化而變化，即同一種截面形狀因子的微通道流動(dòng) Po為一定值。由流體力學(xué)理論可知：對(duì)于常規(guī)不可壓縮充分發(fā)展段流動(dòng)，流動(dòng)阻力系數(shù)與雷諾數(shù)的乘積與雷諾數(shù)無(wú)關(guān)，因此，這一結(jié)果與理論分析相符。可以認(rèn)為Cui等[24]研究中 Po隨雷諾數(shù)的波動(dòng)是實(shí)驗(yàn)測(cè)量誤差等原因造成的。

如圖4所示，Po隨形狀因子的增大而減小，但隨著形狀因子的不斷增大，Po的減小趨勢(shì)逐漸變緩。這表明形狀因子是影響 Po的重要因素，與矩形微通道中的流動(dòng)特征是相似的[6]。在微通道寬度一定的情況下，Po與單位長(zhǎng)度壓降隨形狀因子的變化規(guī)律不同。事實(shí)上，綜合式(1)，(3)和(6)可得：

式中： 為常數(shù)。式(8)表明：當(dāng)寬度與 Re一定時(shí)，微通道單位長(zhǎng)度壓降不僅與阻力系數(shù)相關(guān)，而且與形狀因子相關(guān)；而式(7)則僅取決于阻力系數(shù)，因此，Po與阻力系數(shù)均隨形狀因子的增大而減小，而單位長(zhǎng)度壓降隨形狀因子的變化規(guī)律則有所不同。

圖4 不同Re下微通道形狀因子對(duì)Po數(shù)的影響Fig.4 Effect of micro-channel shape factor on Po number at different Re

3.3 圓角型微通道阻力系數(shù)的經(jīng)驗(yàn)公式

3.3.1 阻力系數(shù)經(jīng)驗(yàn)公式的擬合

截面形狀因子α是影響Po的單一變量。因此，基于圖4中0.01≤Re≤40時(shí)Po隨形狀因子的變化規(guī)律，通過(guò)曲線擬合，得到圓角型微通道阻力系數(shù)的經(jīng)驗(yàn)公式如下：

經(jīng)驗(yàn)公式與數(shù)值模擬結(jié)果之間的誤差不大于0.5%，形狀因子α的適應(yīng)范圍為0.4≤α≤3.2。

3.3.2 阻力系數(shù)經(jīng)驗(yàn)公式對(duì)不同幾何尺寸的微通道適用性

上述研究結(jié)果是在通道截面最大寬度一定、截面高度變化的情況下獲得的，為驗(yàn)證式(9)的普適性，對(duì)截面高度一定(H=50 μm)、截面最大寬度變化的微通道流動(dòng)阻力特性進(jìn)行研究。將基于式(9)計(jì)算得到的不同形狀因子的圓角形微通道 Po預(yù)測(cè)值與對(duì)應(yīng)參數(shù)下的數(shù)值模擬結(jié)果進(jìn)行比較后發(fā)現(xiàn)：當(dāng)0.4≤α≤3.2時(shí)，雷諾數(shù)對(duì)圓角型微通道的 Po沒(méi)有影響，這與前面結(jié)果一致；且式(9)計(jì)算得到的不同形狀因子下的微通道Po預(yù)測(cè)值與數(shù)值模擬結(jié)果之間的誤差小于0.70%。

3.3.3 阻力系數(shù)經(jīng)驗(yàn)公式對(duì)不同黏度流體的適用性

為進(jìn)一步驗(yàn)證阻力系數(shù)經(jīng)驗(yàn)公式的適用性，使用高黏度的離子液(Bmin)BF4對(duì)其進(jìn)行驗(yàn)證，其物理性質(zhì)見(jiàn)表2[25]。表3所示為形狀因子為2.0和不同雷諾數(shù)下基于式(9)的 Po預(yù)測(cè)值與數(shù)值模擬結(jié)果間的相對(duì)誤差。由表3可知：最大誤差小于0.3%。

因此，可以認(rèn)為在所研究的雷諾數(shù)區(qū)域 (0.01≤Re≤40) 和形狀因子范圍(0.4≤α≤3.2)內(nèi)，圓角型微通道阻力系數(shù)經(jīng)驗(yàn)公式(9)具有良好的適用性。

表2 離子液(Bmin)BF4物理性質(zhì)Table 2 Properties of ionic liquid (Bmin)BF4

表3 α=2.0，H=50 μm時(shí)離子液微通道流動(dòng)Po數(shù)相對(duì)誤差分布Table 3 Error distribution of Po of ionic liquid when α=2.0 and H=50 μm

4 結(jié)論

(1) 對(duì)于圓角型微通道，當(dāng)截面寬度一定時(shí)，雷諾數(shù)與截面形狀因子均對(duì)壓降有很大影響，且壓降隨雷諾數(shù)和形狀因子的增大而增大。

(2) 當(dāng)雷諾數(shù)0.01≤Re≤40時(shí)，同一截面形狀因子的圓角形微通道 Po與隨雷諾數(shù)無(wú)關(guān)，因此，截面形狀因子是圓角形微通道 Po的影響單一參變量，且其Po隨截面形狀因子的增大而減小。

(3) 基于圓角型微通道的數(shù)值模擬數(shù)據(jù)，擬合了Poiseuille數(shù)Po隨截面形狀因子變化的阻力系數(shù)α經(jīng)驗(yàn)公式：Po=10.161 1-2.408 5α+0.166 0α2+0.148 09α3-0.197 9α4+0.023 4α5。該式在雷諾數(shù) 0.01≤Re≤40 和截面形狀因子0.4≤α≤3.2范圍內(nèi)均適用。

[1] Jansen H, Gardeniers H, de Boer M, et al. A survey on the reactive ion etching of silicon in microtechnology[J]. Journal of Micromechanics Microengineering, 1996, 6(1): 14-28.

[2] Iliescu C, Tay F E H, MIAO Jian-min. Strategies in deep wet etching of Pyrex glass[J]. Sensors and Actuators A: Physical,2007, 133: 395-400.

[3] Saito Y, Okamoto S, Miki A, et al. Fabrication of microstructure on glass surface using micro-indentation and wet etching process[J]. Applied Surface Science, 2008, 254:7243-7249.

[4] LIU Yong-xun, Ishii K, Tsutsumi T, et al. Ideal rectangular cross-section Si-fin channel double-Gate MOSFETs fabricated using orientation-dependent wet etching[J]. IEEE Electron Device Letters, 2003, 24(7): 484-486.

[5] Iliescu C, MIAO Jian-min, Tay F E H. Optimization of PECVD amorphous silicon process for deep wet etching of Pyrex glass[J].Surface and Coatings Technology, 2005, 192: 43-47.

[6] Shah R K, London A L. Laminar Flow Forced Convection in Ducts[M]. New York: Academic Press, 1978: 157.

[7] Wu P Y, Little W A. Measurement of friction factors for the flow of gases in very fine channels used for microminiaturize Joule–Thomson refrigerators[J]. Cryogenics, 1983, 24(8):273-277.

[8] Harley J, Huang Y, Bau H, et al. Gas flow in microchannels[J].Journal of Fluid Mechanics, 1995, 284: 257-274.

[9] Yang C Y, Chien H T, Lu S R, et al. Friction characteristics of water, R134a and air in small tubes[C]//Celata, G P, ed.Proceedings of the International Conference on Heat Transfer and Transport Phenomena in Microscale. New York: Begell House, 2000: 168-174.

[10] Judy J, Maynes D, Webb B W. Characterization of frictional pressure drop for liquid flows through microchannels[J].International Journal of Heat and Mass Transfer, 2002, 45:3477-3489.

[11] Celata G P, Cumo M, Guglielmi M, et al. Experimental investigation of hydraulic and single phase heat transfer in 0.130 mm capillary tube[J]. Microscale Thermophysical Engineering,2002, 6:85-97.

[12] 馬璨, 袁惠新, 楊振東, 等. 微反應(yīng)器矩形微通道截面高寬比對(duì)流速的影響[J]. 環(huán)境科學(xué)與技術(shù), 2009, 32(10): 64-66.MA Can, YUAN Hui-xin, YANG Zhen-dong, et al. Impact of aspect Ratios on Velocity in Rectangular Micro-channel[J].Environmental Science & Technology, 2009, 32(10): 64-66.

[13] 蔣炳炎, 謝磊, 譚險(xiǎn)峰, 等. 流道截面形狀對(duì)微流體流動(dòng)性能的影響[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2006, 37(5): 964-969.JIANG Bing-yan, XIE Lei, TAN Xian-feng, et al. Influence of section shape of micro channels on microfluid flowability[J].Journal of Central South University:Science and Technology,2006, 37(5): 964-969.

[14] 王瑋, 李志信, 過(guò)增元.粗糙表面對(duì)微尺度流動(dòng)影響的數(shù)值分析[J]. 工熱物理學(xué)報(bào), 2003, 24(1): 85-87.WANG Wei, LI Zhi-xin, GUO Zeng-yuan. Numerical simulation of rough surface effects on microscale fluid flow[J]. Journal of Engineering Thermophysics, 2003, 24(1): 85-87.

[15] 王昊利, 王元, 劉江.平板微通道壁面粗糙度對(duì)流場(chǎng)影響的攝動(dòng)分析[J]. 西安交通大學(xué)學(xué)報(bào), 2005, 39(5): 540-543.WANG Hao-li, WANG Yuan, LIU Jiang. Perturbation analysis of wall roughness effect on flow in micro-channel between two parallel plates[J]. Journal of Xi’an Jiaotong University, 2005,39(5): 540-543.

[16] Hrnjak P, Tu X. Single phase pressure drop in microchannels[J].International Journal of Heat and Fluid Flow, 2007, 28: 1-14.

[17] Yang C Y, Wu J C, Chien H T, et al. Friction characteristics of water, R-134a, and air in small tubes[J]. Microscale Thermophysical Engineering, 2003, 7: 335-348.

[18] 吳嘉峰, 陳永平, 施明恒, 等. 矩形微通道中環(huán)狀冷凝的數(shù)值模擬[J]. 工程熱物理學(xué)報(bào), 2008, 29(11): 1924-1926.WU Jia-feng, CHEN Yong-ping, SHI Ming-heng, et al.Simulation for annular condensation flow in rectangular microchannels[J]. Journal of Engineering Thermophysics, 2008,29(11): 1924-1926.

[19] DUAN Zhi-peng, Muzychka Y S. Slip flow in non-circular microchannels[J]. Microfluidics and Nanofluidics, 2006, DOI 10.1007/s10404-006-0141-4.

[20] QU Wei-lin, Mudawar I. Experimental and numerical study of pressure drop and heat transfer in a single-phase micro-channel heat sink[J]. International Journal of Heat and Mass Transfer,2002, 45: 2549-2565.

[21] Valdés J R, Miana M J, Pelegay J L, et al. Numerical investigation of the influence of roughness on the laminar incompressible fluid flow through annular microchannels[J].International Journal of Heat and Mass Transfer, 2007, 50:1865-1878.

[22] Kandlikar S G, Garimella S, LI Dong-qing, et al. Heat transfer and fluid flow in minichannels and microchannels[M]. Oxford:Elsevier, 2006: 17-21.

[23] Jiang X N, Zhou Z Y, Huang X Y, et al. Laminar flow through microchannels used for microscale cooling systems[C]//Proceedings of the 97 IEEE/CPMT Electronic Packaging Technology Conference. 1997: 119-122.

[24] Cui H H, Silber-Li Z H, Zhu S N. Flow characteristics of liquids in microtubes driven by a high pressure[J]. Physics of Fluids, 2004, 16(5): 1803-1810.

[25] Wilkes J S, Zaworotko M J. Air and water stable 1-ethyl-3-methylimidazolium based ionic liquids[J]. Journal of the Chemical Society, Chemical Communications, 1992(14):965-967.