王琳琳， 李　澤

(1.西安文理學(xué)院 信息工程學(xué)院，西安 710065； 2.西安交通大學(xué) 能源與動(dòng)力工程學(xué)院

熱流科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室，西安 710049； 3.西安工程大學(xué) 理學(xué)院，西安 710048)

?

壁面滑移條件下微尺度通道內(nèi)兩相流數(shù)值模擬

王琳琳1,2， 李澤3

(1.西安文理學(xué)院 信息工程學(xué)院，西安 710065； 2.西安交通大學(xué) 能源與動(dòng)力工程學(xué)院

熱流科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室，西安 710049； 3.西安工程大學(xué) 理學(xué)院，西安 710048)

為研究Y型微通道內(nèi)壁面滑移和連續(xù)相表觀黏度對(duì)通道內(nèi)液滴尺寸和流動(dòng)特性的影響，采用水平集法進(jìn)行數(shù)值模擬.結(jié)果表明：隨著連續(xù)相通道滑移長(zhǎng)度的增加，液滴長(zhǎng)度增長(zhǎng)，通道內(nèi)流動(dòng)更均勻.壁面滑移作用的范圍占通道的46.8%.決定表觀黏性系數(shù)的兩參數(shù)增加，造成負(fù)滑移邊界形成，促進(jìn)液滴脫離，造成液滴長(zhǎng)度縮短.表觀黏度作用范圍占通道的一半以上.對(duì)微流體混合系統(tǒng)的優(yōu)化設(shè)計(jì)提供了參考.

微通道；數(shù)值模擬；水平集法；液滴長(zhǎng)度；邊界滑移；表觀黏度

微型化是當(dāng)今科學(xué)技術(shù)發(fā)展的一個(gè)重要方向.目前，氣泡和液滴在微流體混合器內(nèi)的形成已成為微流動(dòng)動(dòng)力學(xué)研究中的重要問題之一.與常規(guī)尺度通道相比，微流體混合器更利于提高流體間的傳質(zhì)、換熱和反應(yīng)率，具有良好的安全性和可控性[1]，已成功應(yīng)用于微流體乳化[2]、化學(xué)和材料合成、分子進(jìn)化的分割、藥物輸送、藥品檢測(cè)與開發(fā)、生物鑒定、食品工業(yè)和分子篩檢[3]等眾多領(lǐng)域.

常見的微流體混合器分為三類，分別是：T型微通道、匯聚型微通道、同向流動(dòng)微通道[4].毛細(xì)數(shù)表示黏性力和表面張力的比值，毛細(xì)數(shù)小于0.01，微流體混合器的通道壁面約束氣泡或液滴的形成，連續(xù)相的擠壓力是造成其脫離的主要作用力，常見微通道內(nèi)氣泡或液滴的尺寸主要受通道寬度和高度、入口流量比、毛細(xì)數(shù)、壁面條件等因素影響[5-6].目前，對(duì)于其他結(jié)構(gòu)微通道的研究成果相對(duì)較少.宏觀尺度下，通常忽略壁面滑移對(duì)流動(dòng)的影響.微通道內(nèi)，流體在壁面的滑移行為對(duì)流動(dòng)特性產(chǎn)生重要影響，然而對(duì)滑移行為的研究多集中在特定的物理模型，很少將研究滑移問題延伸到微通道內(nèi)兩相流中.

由于微通道尺度微小，測(cè)量微通道內(nèi)氣泡和液滴的形成過程中通道內(nèi)壓降和流場(chǎng)的變化、氣泡和液滴長(zhǎng)度，體積等數(shù)據(jù)對(duì)儀器的精度要求極高，外界極小地?cái)_動(dòng)都可能使實(shí)驗(yàn)結(jié)果發(fā)生顯著變化，測(cè)量設(shè)備價(jià)格高昂.隨著計(jì)算機(jī)運(yùn)行速度和數(shù)值方法精度的提高，數(shù)值模擬已成為研究微流動(dòng)問題的重要手段之一.模擬兩相流界面運(yùn)動(dòng)數(shù)值方法有：標(biāo)記粒子單元法、流體體積法、耗散粒子動(dòng)力學(xué)法、相場(chǎng)法和水平集法等.水平集法是通過引入水平集函數(shù)，兩相界面用該函數(shù)取固定值對(duì)應(yīng)的點(diǎn)集來表示的一種方法，該方法具有捕獲界面精度高、無需重新構(gòu)造界面、能夠準(zhǔn)確描述拓?fù)浣Y(jié)構(gòu)發(fā)生復(fù)雜變化的界面等特點(diǎn)，廣泛應(yīng)用于圖像處理、界面演化、燃燒、流體力學(xué)等領(lǐng)域.

文中基于水平集法，對(duì)Y型結(jié)構(gòu)微通道內(nèi)水-油兩相流動(dòng)進(jìn)行數(shù)值模擬，分析壁面滑移長(zhǎng)度、表觀黏度對(duì)液滴長(zhǎng)度和通道內(nèi)流動(dòng)的影響，為控制微通道內(nèi)液滴尺寸和微流體設(shè)備的優(yōu)化設(shè)計(jì)提供參考.

1　物理模型和控制方程

1.1幾何模型和物性設(shè)定

圖1　Y型微通道幾何模型

文中以Y型微通道為研究對(duì)象，如圖1所示.離散相和連續(xù)相通道之間的夾角，稱為微通道入口角，用α表示，文中幾何模型的微通道入口角是120°.Y型微通道的離散相通道和連續(xù)相通道關(guān)于主通道對(duì)稱，水、油兩相從左側(cè)分別由離散相通道、連續(xù)相通道入口進(jìn)入，在微通道內(nèi)混合后，從右側(cè)出口流出.由此看出，對(duì)流接觸T型微通道是一種特殊的Y型微通道.各通道的寬度都是D1，D1=111μm，離散相和連續(xù)相通道長(zhǎng)度均為3D1，主通道長(zhǎng)度是35D1.水、油兩相密度分別是998kg·m-3、750kg·m-3，黏性系數(shù)分別是0.001Pa·s、0.001 5Pa·s，界面張力系數(shù)是0.047 4N·m-1.

1.2控制方程

微通道內(nèi)水和油互溶性差，壓降小，流體視為不可壓縮的；兩相流體的入口速度較小造成雷諾數(shù)小[7]，通道內(nèi)為層流流動(dòng)；邦德數(shù)小于1，模擬中忽略重力的影響[7].因此，微通道內(nèi)兩相流動(dòng)的連續(xù)性方程和動(dòng)量方程分別是

·U=0

(1)

(2)

式中：U是速度向量，p是壓強(qiáng)，F(xiàn)σ是界面張力，兩相流體密度和黏性系數(shù)

ρ=ρw+(ρo-ρw)φ

(3)

μ=μw+(μo-μw)φ

(4)

下標(biāo)w和o分別代表水和油，兩相流體對(duì)應(yīng)的水平集函數(shù)φ不同，文中水平集函數(shù)采用守恒型，水平集函數(shù)為0和1的區(qū)域分別表示水和油，兩相界面由水平集函數(shù)值為0.5的點(diǎn)組成.水平集函數(shù)滿足方程

(5)

式中：γ是重新初始化參數(shù)，ε是界面厚度.式(2)中的界面張力

Fσ=σκδn

(6)

1.3初始條件和邊界條件

初始時(shí)刻設(shè)定離散相通道充滿水，其他通道充滿油，通道內(nèi)的流體靜止.

微/納米尺度下滑移邊界對(duì)流動(dòng)的影響顯著.微通道內(nèi)的流體主要受壓差的驅(qū)動(dòng)，壁面剪切率低，適合采用Navier滑移模型[8]，即距實(shí)際壁面外側(cè)Ls處為無滑移壁面，Ls稱為滑移長(zhǎng)度.

水、油、壁面接觸時(shí)，分子間引力的作用使相界面上出現(xiàn)浸潤(rùn)現(xiàn)象，油相和壁面間形成接觸角θ，接觸角反映油相對(duì)壁面的吸附能力.接觸角和向量n滿足關(guān)系

n=nwcosθ+twsinθ

(7)

式中：nw、tw分別是壁面單位法向量和切向量.文中接觸角取0°，即油相和通道壁面完全浸潤(rùn).

水油兩相流體分別以u(píng)d=0.017 8m·s-1、uc=0.035 6m·s-1的速度垂直于通道入口流入，通道出口壓強(qiáng)是0Pa，設(shè)置通道壁面是滑移邊界.

文中主要研究Y型微通道內(nèi)兩相流動(dòng)在二維平面的變化特性.采用三角形網(wǎng)格進(jìn)行計(jì)算.通過對(duì)比不同粗細(xì)網(wǎng)格下的計(jì)算結(jié)果，發(fā)現(xiàn)網(wǎng)格數(shù)大于11 000，模擬結(jié)果與網(wǎng)格數(shù)無關(guān)，設(shè)置文中模型網(wǎng)格數(shù)是32 472，時(shí)間步長(zhǎng)取0.000 02s.

2　模擬結(jié)果與分析

2.1滑移長(zhǎng)度對(duì)液滴長(zhǎng)度和連續(xù)相流動(dòng)的影響

速度滑移是指流體在固體表面運(yùn)動(dòng)時(shí)，流體和壁面之間在界面處的切向速度差，此時(shí)壁面為滑移邊界，滑移邊界有助于減小壁面對(duì)流體的拖拽作用.微通道寬度通常小于1mm，壁面對(duì)微流動(dòng)的影響尤其顯著.Choi等通過流變系統(tǒng)測(cè)得水的滑移長(zhǎng)度為2～20μm，體積分?jǐn)?shù)30%甘油的滑移長(zhǎng)度為5～50μm.劉趙淼等數(shù)值模擬研究了油膜厚度為50μm、滑移長(zhǎng)度為0.5～10μm的微流動(dòng).為進(jìn)一步研究Y型微通道內(nèi)滑移邊界對(duì)液滴長(zhǎng)度和流動(dòng)的影響規(guī)律，取離散相通道壁面滑移長(zhǎng)度為1μm，連續(xù)相和主通道的滑移長(zhǎng)度在1～20μm之間變化.隨著滑移長(zhǎng)度的增加，液滴長(zhǎng)度L逐漸增加，但滑移長(zhǎng)度由14μm繼續(xù)增加，液滴長(zhǎng)度增長(zhǎng)減緩，如圖2所示.

液滴后端油相的擠壓力對(duì)液滴表面起破壞作用，擠壓力增大促進(jìn)液滴脫離.為研究滑移長(zhǎng)度對(duì)油相擠壓作用的影響，在連續(xù)相通道內(nèi)取點(diǎn)A1，見圖1，取滑移長(zhǎng)度分別為4、8、12、16和20μm，液滴形成過程中點(diǎn)A1壓強(qiáng)的變化如圖3所示.圖中可見，滑移長(zhǎng)度增加，油相壓強(qiáng)減?。换崎L(zhǎng)度由14μm繼續(xù)增大，壓強(qiáng)減小程度變?nèi)?這些結(jié)果表明，滑移長(zhǎng)度增加造成液滴后端受到油相的擠壓作用減弱，破壞力的減小使液滴脫離時(shí)間延長(zhǎng)，從而造成液滴長(zhǎng)度隨壁面滑移程度的增加而增加；然而，隨著滑移長(zhǎng)度增加，油相對(duì)液滴后端擠壓力減小程度變?nèi)酰斐梢旱伍L(zhǎng)度的增長(zhǎng)程度減弱，由此解釋了圖2中模擬結(jié)果的原因.

圖2　液滴長(zhǎng)度隨壁面滑移長(zhǎng)度的變化

圖3　連續(xù)相通道內(nèi)壓強(qiáng)隨滑移長(zhǎng)度的變化

為研究滑移長(zhǎng)度對(duì)通道內(nèi)速度分布的影響，在連續(xù)相通道內(nèi)建立坐標(biāo)系，如圖4(a)，x1軸正方向?yàn)橛拖嗳肟诜较?，y1軸垂直于該方向.連續(xù)相通道內(nèi)油相的流動(dòng)是泊肅葉流動(dòng)，沿x1軸正方向的流速uo滿足

(8)

結(jié)合滑移壁面條件對(duì)式(8)積分，得出油相的流速

(9)

圖4(b)顯示了不同滑移長(zhǎng)度下，通過數(shù)值模擬和理論值式(9)得到的連續(xù)相通道內(nèi)油相沿徑向的速度分布.圖中可見數(shù)值模擬和理論公式得到的油相速度幾乎相同，速度沿連續(xù)相通道軸線呈拋物線對(duì)稱分布，各速度曲線交于公共點(diǎn)，交點(diǎn)分別在y1=0.026mm和y1=0.083mm，說明兩交點(diǎn)距壁面占通道的46.8%左右，得到連續(xù)相通道內(nèi)壁面滑移作用的范圍.隨滑移長(zhǎng)度增加，油相速度在通道中心線上取得的最大值相對(duì)于油相入口速度由1.348減至1.24，而壁面上的滑移速度相對(duì)于油相入口速度由0.301增至0.519，這些結(jié)果表明，滑移長(zhǎng)度的增加促進(jìn)了通道內(nèi)連續(xù)相均勻流動(dòng).

圖4　連續(xù)相通道內(nèi)建立的x1-y1坐標(biāo)系和油相速度分布隨滑移長(zhǎng)度的變化

由滑移壁面條件知，連續(xù)相通道壁面對(duì)油相施加的黏性切應(yīng)力

(10)

圖5　壁面切應(yīng)力隨滑移長(zhǎng)度的變化

式中：us是壁面的滑移速度.不同滑移長(zhǎng)度下，由式(10)得到壁面相對(duì)切應(yīng)力如圖5所示.圖中可見，滑移長(zhǎng)度增加，油相受到壁面黏性切應(yīng)力減小，減小程度隨滑移長(zhǎng)度的增加而減弱，滑移長(zhǎng)度由14 μm繼續(xù)增加，該切應(yīng)力隨滑移長(zhǎng)度由二次減小變?yōu)榫€性減小.由此表明，滑移長(zhǎng)度增加，造成油相受到壁面的黏性阻力減弱，即滑移壁面具有減阻的效果，但滑移壁面減阻作用會(huì)隨著滑移長(zhǎng)度的增加而減弱.

2.2固壁分子作用對(duì)液滴長(zhǎng)度的影響

宏觀尺度下流體黏性稱為常規(guī)黏性，由分子間的吸引力決定，其大小與分子間吸引力成正比.微尺度流動(dòng)中，液體黏性除常規(guī)黏性外，還有一部分是由壁面和近壁面處液體分子間相互吸引而產(chǎn)生的附加黏性，受壁面影響的液體黏性系數(shù)為表觀黏性系數(shù)μa

μa=μ+Δμ

(11)

式中：μ為宏觀尺度下流體的黏性系數(shù)，Δμ為附加黏性系數(shù).

附加黏性系數(shù)隨流體和壁面間距減小而增加.對(duì)于特征尺度較大的流動(dòng)通常忽略壁面對(duì)流體分子吸引產(chǎn)生的黏性.微尺度的流動(dòng)中不能忽略這種壁面的吸引作用，在二維微流動(dòng)中，考慮到通道壁面對(duì)流體的共同作用，將附加黏性系數(shù)定義為

(12)

式中無量綱數(shù)ξ和n由流體和固壁確定，d1是流體質(zhì)點(diǎn)到通道軸線的距離，為了防止壁面上流體黏度無限大，取a的值為1.08.可以看出，附加黏性系數(shù)受n和ξ影響.分別改變兩參數(shù)值，以詳細(xì)研究表觀黏性對(duì)微通道內(nèi)液滴尺寸和流動(dòng)的影響規(guī)律.

2.2.1n對(duì)液滴長(zhǎng)度和連續(xù)相流動(dòng)的影響

為研究參數(shù)n對(duì)液滴長(zhǎng)度的影響，連續(xù)相通道滑移長(zhǎng)度是8 μm的通道內(nèi)，固定無量綱數(shù)ξ為0.1，n分別取1、2、3，得到長(zhǎng)度分別為3.02D1、3.01D1、2.94D1的液滴.由此可見，參數(shù)n增加造成液滴長(zhǎng)度縮短.

圖6　連續(xù)相通道內(nèi)速度分布隨參數(shù)n的變化

對(duì)于無附加黏度和不同n對(duì)應(yīng)的附加黏度，連續(xù)相通道y1軸上油相的速度分布變化如圖6所示.圖中可見，無附加黏度，速度分布呈拋物線狀；隨著n增加，壁面上的速度由0.218uc減小至0.008uc，連續(xù)相通道內(nèi)油相速度最大值明顯增加，由1.449uc增至2.037uc，速度分布不均勻程度增加，油相流動(dòng)越來越集中在連續(xù)相通道的中間區(qū)域.圖中的速度曲線基本在相同位置相交，交點(diǎn)分別位于y1=0.03 mm和y1=0.081 mm，交點(diǎn)至壁面的距離占通道的54.1%左右，得到了連續(xù)相通道內(nèi)受參數(shù)n影響的油相流體表觀黏性作用范圍.

圖中顯示n大于1，油相速度分布曲線上出現(xiàn)拐點(diǎn)，兩拐點(diǎn)間的區(qū)域?yàn)橛拖嘀髁鲄^(qū)，油相在主流區(qū)速度分布呈拋物線狀，在拐點(diǎn)處做主流區(qū)速度分布曲線的切線，該切線和縱軸負(fù)半軸相交，表明此時(shí)連續(xù)相通道壁面上的滑移速度為負(fù)值，定義此邊界為負(fù)滑移邊界，它在模擬宏觀和納米尺度的流動(dòng)中有較多應(yīng)用.圖中可見，n增大，拐點(diǎn)逐漸遠(yuǎn)離通道壁面，造成油相在壁面上負(fù)滑移的程度增大，油相的主流區(qū)變窄.速度分布曲線上拐點(diǎn)的出現(xiàn)，導(dǎo)致油相的流動(dòng)容易失穩(wěn).n越大，越容易使流動(dòng)失穩(wěn)，造成n取3時(shí)，形成的液滴長(zhǎng)度最短.

2.2.2ξ對(duì)液滴長(zhǎng)度和連續(xù)相流動(dòng)的影響

為考察ξ對(duì)液滴長(zhǎng)度的影響，固定n為2，連續(xù)相通道壁面滑移長(zhǎng)度為8 μm，參數(shù)ξ取不同值得到不同長(zhǎng)度的液滴，如圖7所示.圖中可見，隨著ξ由0增至0.1，液滴長(zhǎng)度由3.09D1減小至3.01D1；ξ由0.1增加到0.5，液滴長(zhǎng)度僅由3.01D1減小至2.99D1；ξ繼續(xù)增加，液滴長(zhǎng)度幾乎不變.

圖7　液滴長(zhǎng)度隨參數(shù)ξ的變化

圖8　連續(xù)相通道內(nèi)速度分布隨參數(shù)ξ的變化

改變?chǔ)卧斐蛇B續(xù)相通道內(nèi)y1軸上油相速度分布發(fā)生變化，如圖8所示.圖中可見，速度分布曲線基本在y1=0.032 mm和y1=0.079 mm處相交，表明連續(xù)相通道內(nèi)受ξ影響的油相流體表觀黏性作用范圍是57.7%，比參數(shù)n的影響范圍略大.隨著ξ增加，油相在壁面上的速度由0.301uc減小至0.028 4uc，通道中心線上的最大速度由1.35uc增至1.972uc，但速度分布曲線的變化逐漸減弱，當(dāng)ξ由0.5增至0.6，速度分布幾乎不變.ξ大于零，速度分布曲線出現(xiàn)拐點(diǎn)，隨著ξ增大，拐點(diǎn)遠(yuǎn)離通道壁面，造成油相在壁面上負(fù)滑移程度增長(zhǎng)，油相的主流區(qū)變窄，但ξ由0.5增至0.6，油相的負(fù)滑移程度幾乎不變，由此造成ξ由0.5繼續(xù)增加，液滴長(zhǎng)度基本相同.

3　結(jié)語

基于水平集法，對(duì)Y型微通道內(nèi)水-油兩相流動(dòng)進(jìn)行數(shù)值模擬，通過改變連續(xù)相通道的滑移長(zhǎng)度和油相表觀黏度，得到下面的結(jié)果：

1)滑移長(zhǎng)度增加，促進(jìn)了液滴長(zhǎng)度和壁面滑移速度的增加，使連續(xù)相通道內(nèi)的速度分布更加均勻.連續(xù)相通道內(nèi)壁面滑移作用范圍占通道的46.8%.

2)增加影響油相表觀黏度的參數(shù)n或ξ，液滴長(zhǎng)度均縮短.參數(shù)n大于1或ξ大于零，連續(xù)相通道出現(xiàn)負(fù)滑移邊界，通道內(nèi)流動(dòng)失穩(wěn)，促進(jìn)液滴脫離.

3)連續(xù)相通道內(nèi)受參數(shù)n影響的油相流體表觀黏性作用范圍占連續(xù)相通道的54.1%，受參數(shù)ξ影響的作用范圍略大，占57.7%.

[1]WANGWei,XIERui,JUXiao-jie,etal.Controllablemicrofluidicproductionofmulticomponentmultipleemulsions[J].LabonaChip,2011,11(9):1587-1592.

[2]JANKOWSKIP,OGOCZYKD,DERZSIL,etal.Hydrophilicpolycarbonatechipsforgenerationofoil-in-water(O/W)andwater-in-oil-in-water(W/O/W)emulsions[J].MicrofluidicsandNanofluidics,2013,14(5): 767-774.

[3]PAWELJ,DOMINIKAO,LADISLAVD,etal.Hydrophilicpolycarbonatechipsforgenerationofoil-in-water(O/W)andwater-in-oil-in-water(W/O/W)emulsions[J].MicrofluidNanofluid,2013,14(1):597-604.

[4]ZHAOChun-xia,ANTONPJ.Two-phasemicrofluidicflows[J].ChemicalEngineeringScience,2011,66:1394-1411.

[5]VANSTEIJNV,KLEIJNCR,KREUTZERMT.PredictivemodelforthesizeofbubblesanddropletscreatedinmicrofluidicT-junctions[J].LabChip,2010,10:2513-2518.

[6]WANGXiao-da,ZHUChun-ying,WUYi-ning,etal.DynamicsofbubblebreakupwithpartlyobstructioninamicrofluidicT-junction[J].ChemicalEngineeringScience,2015,132:128-138.

[7]袁希鋼,宋文琦.T型結(jié)構(gòu)微通道氣液兩相流型的數(shù)值模擬[J].天津大學(xué)學(xué)報(bào),2012,45(9):763-769.

[8]劉趙淼,王國斌,申峰.基于Navier滑移的油膜縫隙微流動(dòng)特性數(shù)值分析[J].機(jī)械工程學(xué)報(bào)，2011,47(21):104-110.

[責(zé)任編輯馬云彤]

Numerical Simulation of the Two-phase Flow in a Micro-scale Channelunder the Boundary Slip Condition

WANG Lin-lin1,2, LI Ze3

(1.SchoolofInformationEngineering,Xi’anUniversity,Xi’an710065,China; 2.KeyLaboratoryofThermalFluidScienceandEngineeringofMOE,SchoolofEnergyandPowerEngineering,Xi’anJiaotongUniversity,Xi’an710049,China; 3.SchoolofScience,Xi’anPolytechnicUniversity,Xi’an710048,China)

InordertoinvestigatetheinfluenceofthewallslipandapparentviscosityofthecontinuousphaseonthedropletsizeandflowcharacteristicsinaY-junctionmicro-channel,thelevelsetmethodisusedfornumericalsimulation.Theresultsshowthatwhenthesliplengthincreasesinthecontinuousphasechannel,thedropletlengthincreases,andthevelocitydistributionbecomesevenly.Therangeoftheboundaryslipeffectaccountsfor46.8%ofthechannel.Theapparentviscosityisdeterminedbythetwoparameters.Thenegativeboundaryslipisformed,thedropletdetachmentisaccelerated,anddropletlengthdecreasesforincreasingthetwoparameters.Therangeoftheapparentviscosityeffectaccountsformorethanhalfofthechannel.Itprovidesareferencefortheoptimaldesignofthemicrofluidmixingsystem.

micro-channel;numericalsimulation;levelsetmethod;dropletlength;boundaryslip;apparentviscosity

1008-5564(2016)03-0028-06

2016-01-04

國家自然科學(xué)基金資助項(xiàng)目(51076126);西安市科技計(jì)劃項(xiàng)目(CXY1443WL19)

王琳琳(1981—)，女，河南長(zhǎng)垣人，西安文理學(xué)院信息工程學(xué)院講師，博士，主要從事微通道內(nèi)兩相流動(dòng)研究.

TQ021.1

A