類文丘里十字型流道中的液滴流態(tài)分析

夏侯彬,蔣志韜,陳忠利,錢麗娟,2

(1.中國計量大學 機電工程學院,浙江 杭州 310018;2.浙江省智能制造質(zhì)量大數(shù)據(jù)溯源與應(yīng)用重點實驗室,浙江 杭州 310018)

微液滴以其尺寸小,可控性強,比表面積大等優(yōu)點被廣泛應(yīng)用于生物醫(yī)學[1]、藥物篩選和輸送[2]、食品工業(yè)[3]和化妝品制造[4]等行業(yè)。通常情況下,利用微流控技術(shù)生成微液滴主要依賴于多相流的流體流動,但同時也與所選擇的微流道幾何形狀密切相關(guān)。這是由于改變微流道的幾何形狀可以有效地控制其內(nèi)部多相流的流體動力學和界面演化過程,從而更好地制備所需特性的微液滴。

科研人員采用實驗和數(shù)值模擬方法,就十字型流道結(jié)構(gòu)變化對液滴生成的影響做了一系列研究。Liu和Zhang利用格子玻爾茲曼方法(lattice Boltzmann method, LBM)對十字型微流道進行了三維數(shù)值模擬,發(fā)現(xiàn)當微流道的寬深比以及兩相進口的寬度比越大時,生成的液滴尺寸也就越大[5]。Gulati等通過實驗探究得出十字交叉處的圓形結(jié)構(gòu)會對液滴生成產(chǎn)生阻礙作用[6]。Zhang等利用水平集(level set method, LS)方法,比較了十字交叉型出口處4種不同頸縮結(jié)構(gòu)下液滴生成的尺寸和周期,發(fā)現(xiàn)當出口頸縮結(jié)構(gòu)為梯形時生成的液滴尺寸最小,頻率最高[7]。Ngo和Iqbal等采用流體體積法(volume of fluid method,VOF)對十字流道進口夾角進行了研究,觀測得出液滴尺寸在夾角為90°時最大,而隨著進口夾角偏離90°,液滴尺寸會逐漸減小[8-9]。

綜上所述,十字型流道結(jié)構(gòu)變化的研究多集中在結(jié)構(gòu)變化對于液滴尺寸和生成周期的影響。然而,對于流道結(jié)構(gòu)變化下液滴流態(tài)的分析以及液滴生成時壓力積聚現(xiàn)象并未有太多提及。查閱相關(guān)文獻發(fā)現(xiàn),文丘里結(jié)構(gòu)可以有效緩解流道內(nèi)的壓力積聚問題,已有前人將文丘里結(jié)構(gòu)加入T型流道中,來減緩液滴生成時的壓力積聚問題[10]。因此,本文采用VOF方法,對類文丘里十字型流道的液滴生成流態(tài)進行分析,以便更好地理解液滴生成的機制,從而為后續(xù)實現(xiàn)液滴精確控制提供一定指導。

1 數(shù)值模型

1.1 控制方程

VOF界面追蹤方法[11]利用體積函數(shù)α表示某一項流體在一個單元網(wǎng)格中所占體積分數(shù)。本文以分散相體積為主要研究對象,每個網(wǎng)格單元中兩相流體的體積分數(shù)關(guān)系如下:

αd+αc=1。

(1)

式(1)中下標d和c分別代表分散相流體和連續(xù)相流體。

通過求解體積函數(shù)的輸運方程來追蹤兩相界面,如式(2):

(2)

對于單元網(wǎng)格中的αd來說,若αd=1,則表示單元網(wǎng)格中全為分散相,若αd=0,則代表單元網(wǎng)格中全為連續(xù)相,若0<αd<1,則在單元網(wǎng)格中既有分散相又有連續(xù)相,即存在兩相的交界面。

假設(shè)兩相流體均為不可壓縮的牛頓流體,其滿足的連續(xù)性方程和動量方程如式(3)和式(4):

▽·U=0,

(3)

(4)

式(3)和(4)中:U是流體速度矢量,m/s;μ是流體動力粘度,Pa·s;ρ是流體密度,kg/m3;P是壓力,Pa;g是重力加速度,m/s2;Fs是兩相間的界面張力,N/m3。μ和ρ是流體的體積加權(quán)平均數(shù),可以通過兩相流體的體積分數(shù)表示,如下:

μ=μdαd+(1-αd)μc,

(5)

ρ=ρdαd+(1-αd)ρc。

(6)

Fs表示界面張力,可以通過式(7)的連續(xù)界面張力模型(continuum surface force,CSF)[12]計算得出

(7)

式(7)中,σ為表面張力系數(shù),N/m;κ為界面曲率,可由式(8)計算得到

(8)

(9)

當相界面與通道內(nèi)壁接觸時,需要考慮壁面的粘附效應(yīng)[13],它通常由壁面接觸角決定,且接觸角用于計算壁面附近的界面曲率,如式(10):

(10)

1.2 數(shù)值模擬模型

本文所用的十字型微流道的二維數(shù)值模擬模型如圖1,兩相入口寬度和長度均相等,通道寬度設(shè)置為Wc=Wd=100 μm,通道長度設(shè)置為Lc=Ld=100 μm。將交匯處改為類文丘里結(jié)構(gòu),其中縮進長度設(shè)置為Wv=Hv=25 μm,縮進角度設(shè)置為θ=45°。為了確保兩相能夠在主流道中充分流動,將主流道長度設(shè)置為Lmain=1 500 μm。分散相以流量Qd注入分散相入口,連續(xù)相以流量Qc注入對稱的連續(xù)相入口。為簡化計算成本,將連續(xù)相和分散相的物性參數(shù)設(shè)置為ρd=ρc=1 000 kg/m3、μd=0.001 Pa·s、μc=0.005 Pa·s,兩相間界面張力系數(shù)設(shè)置為σ=0.01 N/m,接觸角設(shè)置為θw=150°。

圖1 十字型流道數(shù)值模擬結(jié)構(gòu)示意圖Figure 1 Schematic diagram of numerical simulation flow-focusing device

1.3 邊界條件與數(shù)值方法

為了簡化計算模型,對計算模型進行如下假設(shè):兩相流動均為不可壓縮流動;與外界無熱量交換;分散相與連續(xù)相在流動過程中無化學反應(yīng)發(fā)生;兩相流體的物性參數(shù)為常數(shù)。由于微流道流動中的雷諾數(shù)小于0.1,故選擇層流模型和基于壓力基的非穩(wěn)態(tài)求解器。邊界條件設(shè)置如下:兩相入口都設(shè)置為恒定的速度入口,出口設(shè)置為壓力出口,出口壓力Pout=0 Pa。

使用VOF框架中的有限體積法求解控制方程以及壓力隱式分裂算子(PISO)求解連續(xù)性方程和動量方程。壓力方程采用PRESTIO!算法離散,而動量方程則通過二階迎風格式離散。相界面重構(gòu)采用分段線性界面法[14]。在全局庫朗數(shù)小于0.2的標準下,采用變時間步長法確保計算的收斂和穩(wěn)定。松弛因子設(shè)定如下:壓力0.2,密度0.3,體積力0.3,動量0.2。

1.4 網(wǎng)格獨立性驗證

為了尋找合適的網(wǎng)格尺寸,本章采用結(jié)構(gòu)化四邊形網(wǎng)格對微流道的二維幾何模型進行劃分,網(wǎng)格劃分如圖2。選取三種不同網(wǎng)格數(shù)進行網(wǎng)格獨立性驗證,即32 400、41 560和64 750。定義兩個形貌參數(shù)來進行獨立性驗證,分別是分散相的液線長度l和分散相的頸部寬度δ(見圖1下方插圖)。其中,頸部寬度在膨脹階段定義為縱向最大寬度,在頸縮階段定義為頸部的最小寬度。

圖2 二維十字型流道網(wǎng)格示意圖Figure 2 Schematic diagram of the grid of two-dimensional flow-focusing microchannels

圖3(a)為膨脹和頸縮階段下某一時刻的液滴形貌圖。圖3(b)為分散相液線長度隨時間的演化過程。從圖3可以看出,隨著網(wǎng)格數(shù)的增加,數(shù)值模擬結(jié)果逐漸收斂。當網(wǎng)格數(shù)為41 560和64 750時,數(shù)值模擬結(jié)果趨于一致,考慮到計算和時間成本,故選用網(wǎng)格數(shù)為41 560的網(wǎng)格尺寸。

圖3 不同網(wǎng)格數(shù)N=32 400、41 560和64 750下網(wǎng)格獨立性檢驗結(jié)果Figure 3 Results of grid independence tests for different grid numbers N=32 400, 41 560 and 64 750

1.5 模型驗證

為了驗證模型的正確性,選取Soroor[15]所做實驗進行對應(yīng)的數(shù)值模擬。實驗采用的十字型微流道的主通道寬度為125 μm交匯處出口的頸縮寬度為50 μm,長度為125 μm。選取甘油-水溶液(90%的水,10%的甘油)作為分散相(ρd=1 021.90 kg/m3,μd=1.16 mPa·s),礦物油作為連續(xù)相(ρc=830.66 kg/m3,μc=38.5 mPa·s),界面張力系數(shù)為σ=37.1 mN/m。數(shù)值模擬驗證的實驗條件為Qd=8 μL/min,Qc=7.62 μL/min。

數(shù)值模擬所得結(jié)果和實驗對比情況如圖4,左邊為實驗圖,圖中紅色線條為采用LS方法所得的對照結(jié)果,右圖為采用VOF方法所得的驗證結(jié)果,從圖4可以看出,盡管數(shù)值模擬中某些時刻液滴的輪廓存在差異,但整體曲率輪廓和文獻中實驗所得的曲率輪廓幾乎吻合,故可以保證數(shù)值模型的正確性。

圖4 數(shù)值模擬和實驗所得的液滴生成過程的形貌演化對比圖Figure 4 Morphological evolution of droplet generation process obtained from numerical simulation and experiment

2 流態(tài)分析

2.1 擠壓態(tài)分析

數(shù)值模擬結(jié)果表明,擠壓態(tài)出現(xiàn)在Ca較小且Q較大處。圖5展示了擠壓態(tài)在一個生成周期內(nèi)相界面形貌、流場和壓力場的非穩(wěn)態(tài)演化。根據(jù)液滴生成過程中相界面的演化過程,將擠壓態(tài)液滴生成過程分為4個階段:膨脹階段(expanding),擠壓階段(squeezing),頸縮階段(necking)和斷裂階段(pinch-off)。

膨脹階段(見圖5,t=0～3.1 ms):從流線流動方向發(fā)現(xiàn)連續(xù)相大部分沿著交匯處的空隙流出,表明分散相受連續(xù)相的擠壓作用較弱。隨著交匯處分散相體積不斷增大,連續(xù)相流出空間減小,造成交匯處壓力積聚和升高。

擠壓階段(見圖5,t=3.1～6.5 ms):連續(xù)相的擠壓作用增強,分散相的徑向?qū)挾乳_始減小。分散相頭部進入主流道,形成柱塞狀。從流場變化可知,連續(xù)相的流線從交匯處的出口轉(zhuǎn)向進口,流線的聚集程度減緩,使得交匯處的壓力減小。

頸縮階段(見圖5,t=6.5～10.0 ms):分散相的頸部內(nèi)凹變形并逐漸拉長變細,以抵抗連續(xù)相的擠壓作用。連續(xù)相的流線垂直于分散性的流線,促進頸部的凹陷,使得頸部處的壓力增大。

斷裂階段(見圖5,t=10.0～11.5 ms):連續(xù)相流線的垂直分布表面連續(xù)相對頸部的擠壓作用仍然存在,使得壓力繼續(xù)增大,促進頸部快速斷裂形成液滴。在表明張力作用下,液滴尾部與斷裂后的分散相前端后退回縮。

為了更好的理解液滴生成的機理,采用1.4節(jié)中定義的分散相的液線長度l和頸部寬度δ以及它們各自的變化率(Δl/Δt,Δδ/Δt)定量描述液滴生成過程中相界面的形貌演化。此外,將分散相頭部前端界面內(nèi)外的拉普拉斯壓差(Δp=pi-po),見圖5(a),t=0 ms,隨時間的變化作為考察對象。從圖5可以看出,在膨脹階段,分散相頭部在兩相交匯處體積持續(xù)擴大,l和δ逐漸增大,見圖6(a)(c),t=0～3.1 ms。該階段內(nèi)來自側(cè)通道的連續(xù)相對交匯處的分散相的擠壓作用增強,進而抑制了分散相頭部在y方向上的膨脹,因此Δδ/Δt減小,見圖6(d),t=0～3.1 ms。而這又促進了分散相頭部向主流道下游的移動,使得Δl/Δt增大,見圖6(b),t=0～3.1 ms。而后,當分散相頭部到達兩相交匯處出口后,由于流出空間的減小,分散相在y方向?qū)B續(xù)相的擠壓作用的抵抗增強,故Δl/Δt緩慢減小,見圖6(b),t=0～3.1 ms,Δδ/Δt緩慢增大,見圖6(d),t=0～3.1 ms。在膨脹階段初期,分散相頭部在慣性力的作用下不斷膨脹,前端界面的曲率半徑增大,Δp隨之減小,見圖6(e),t=0～3.1 ms,液滴穩(wěn)定性增加,進而液滴前端界面趨向圓形。而在膨脹階段后期,連續(xù)相的擠壓作用加上流出空間的減小,使得分散相頭部前端界面不斷前凸,界面曲率半徑減小,進而Δp增大,見圖6(e),t=0～3.1 ms,液滴穩(wěn)定性減弱,液滴前端界面偏離圓形。

在擠壓階段,隨著連續(xù)相擠壓作用的持續(xù)增強,交匯處內(nèi)凸起的分散相受到擠壓,δ開始縮小,Δδ/Δt變?yōu)樨撝登移浣^對值隨著時間持續(xù)增大,表明δ的減小速率不斷上升,見圖6(c)(d),t=3.1～6.5 ms。同時,l與Δl/Δt隨時間明顯增大,見圖6(a)(b),t=3.1～6.5 ms,說明連續(xù)相的擠壓作用還使得分散相頭部向下游延伸。隨著分散相從兩相交匯處被擠入主通道下游的有限空間內(nèi),頭部前端的曲率半徑快速減小,Δp迅速增大,見圖6(e),t=3.1～6.5 ms,液滴穩(wěn)定性進一步減弱,液滴前端界面由圓形變?yōu)榧忓N形。

從擠壓階段過渡到頸縮階段時,兩相交匯處被分散相完全堵塞,連續(xù)相的擠壓作用達到最大,頸部開始內(nèi)凹,Δδ/Δt達到最小,見圖6(c),t=6.5～10.0 ms,并趨于穩(wěn)定。而且,隨著頸部收縮程度和分散相頭部體積的增大,頸縮過程對l的促進作用大幅降低。因此,該階段內(nèi),l隨恒定流量的分散相的注入近似表現(xiàn)出線性增長,Δl/Δt小幅度的減小最后趨于穩(wěn)定,見圖6(a)(b),t=6.5～10.0 ms。隨著主流道分散相頭部的逐漸穩(wěn)定,界面半徑趨向穩(wěn)定,Δp緩慢減小,見圖6(e),t=6.5～10.0 ms,液滴前端界面又回到圓形。

在斷裂階段,頸部在Rayleigh-Plateau不穩(wěn)定性[16]的誘導下斷裂,δ減小到0,由于分散相對出流口的堵塞,Δδ/Δt并無較大變化,見圖6(c)和(d)。由于斷裂處頸部的體積很小,液滴頸部的夾斷對分散相液線向下游移動的影響極其有限,因此,在該階段內(nèi),分散相的液線長度延續(xù)了頸縮階段結(jié)束時的變化率,即Δl/Δt和Δδ/Δt保持不變,見圖6(b)。此外分散相頭部前端的界面形貌已經(jīng)達到穩(wěn)定,故Δp基本不變,見圖6(e),液滴前端界面保持圓形不變。

2.2 滴落態(tài)分析

在擠壓態(tài)的基礎(chǔ)上,減小流量比和毛細數(shù),液滴生成流態(tài)轉(zhuǎn)變?yōu)榈温鋺B(tài)。滴落態(tài)下液滴相界面形貌、流場和壓力場的非穩(wěn)態(tài)演化如圖7。與擠壓態(tài)相同,滴落態(tài)中液滴生成也分為4個階段。

圖7 滴落態(tài)液滴流體動力學行為的演化過程Figure 7 Evolution of the fluid dynamics behavior of dripping regime

膨脹階段(見圖7,t=0～13.25 ms):分散相的流量減小,交匯處分散相體積增大速率變慢,流線分布呈現(xiàn)出中間密兩端疏的狀態(tài),兩相交匯處的壓力緩慢上升。

擠壓階段(見圖7,t=13.25～18.50 ms):連續(xù)相粘性剪切作用下交匯處分散相開始收縮,主流道中分散相頭部為圓形。連續(xù)相的流線從交匯處的出口轉(zhuǎn)向進口,聚集程度減緩,交匯處的壓力略微減小。

頸縮階段(見圖7,t=18.50～21.50 ms):頸部受到連續(xù)相粘性剪切作用開始內(nèi)凹變薄,主流道內(nèi)分散相體積緩慢增大。連續(xù)相流線垂直于頸部,流線匯聚使得壓強升高。

斷裂階段(見圖7,t=21.50～23.75 ms):在連續(xù)相粘性剪切作用下,流線匯聚位置相較于擠壓態(tài)更靠近分散相入口,進而交匯處壓強繼續(xù)增大,促進頸部斷裂形成液滴。

滴落態(tài)液滴界面演化特征參數(shù)如圖8。在膨脹階段,連續(xù)相的粘性剪切作用作為主導[17],l和δ平穩(wěn)增加,Δl/Δt和Δδ/Δt保持不變,見圖8(a)～(d),t=0～13.25 ms。此外,分散相頭部前端的曲率半徑緩慢增加,Δp緩慢減小,見圖8(e),t=0～13.25 ms。

圖8 滴落流態(tài)下液滴生成過程中特征參數(shù)演化Figure 8 Evolution of characteristic parameters during droplet generation in dripping regime

進入擠壓階段,連續(xù)相的粘性剪切作用持續(xù),l和δ繼續(xù)保持穩(wěn)定增長,Δl/Δt和Δδ/Δt保持不變,見圖8(a)～(d),t=13.25～18.50 ms。在擠壓階段后期,由于分散相頭部體積的增大,使得連續(xù)相的粘性剪切與拖曳作用持續(xù)增強,進而引起頸部變細速率和頭部伸長速率的加快,因此Δl/Δt和Δδ/Δt的絕對值隨時間增大,見圖8(b)和(d),t=13.25～18.50 ms。由于頭部體積的增大,曲率半徑開始減小,故Δp增大,見圖8(e),t=13.25～18.50 ms。在頸縮階段,分散相頭部尺寸趨于穩(wěn)定使得來自連續(xù)相的粘性剪切和拖曳作用的增強速率也隨之降低并趨于穩(wěn)定,進而引起Δl/Δt和Δδ/Δt的絕對值的增大逐漸變緩并趨于穩(wěn)定,見圖8(b)和(d),t=18.50～21.50 ms。分散相頭部尺寸的增大,使得前端的曲率半徑繼續(xù)增大,進而Δp減小,見圖8(e),t=18.50～21.50 ms。

斷裂階段,Rayleigh-Plateau不穩(wěn)定性引起液滴頸部斷裂,δ變?yōu)?,Δδ/Δt的絕對值急劇增大,見圖8(c)和(d)。不同于擠壓態(tài),盡管液滴頭部體積達到穩(wěn)定,但此時分散相頭部前端在界面張力的作用下逐漸變?yōu)榍蛐?曲率半徑不斷增大,Δp繼續(xù)減小,見圖8(e)。

2.3 射流態(tài)流態(tài)分析

在擠壓態(tài)的基礎(chǔ)上,保持毛細數(shù)不變繼續(xù)增大流量比,液滴流態(tài)轉(zhuǎn)變?yōu)樯淞鲬B(tài)。相較于擠壓態(tài)和滴落態(tài),射流態(tài)下并無液滴生成,而是形成一條均勻穩(wěn)定的液線,隨著兩側(cè)通道的連續(xù)相一起沿主流道流出,如圖9。在射流態(tài)下,分散相的慣性力主導分散相液線的流動,且足夠大的分散相慣性力可以抑制在兩相交界處由界面張力誘發(fā)的Rayleigh-Plateau不穩(wěn)定性,因此分散相與連續(xù)相界面處不會發(fā)生波動和破裂。

圖9 射流態(tài)流體動力學行為的演化過程Figure 9 Evolution of the fluid dynamics behavior of threading regime

3 流態(tài)圖分析

如上所述,毛細數(shù)和流量比的變化改變了液滴生成時不同力之間的競爭關(guān)系。因此,基于毛細數(shù)和流量比,繪制了流態(tài)分布圖,如圖10。當毛細數(shù)和流量比都較小時,連續(xù)相的粘性力和分散相的慣性力也很小,分散相在界面張力的主導下在兩相交匯處經(jīng)歷膨脹和擠壓過程,即出現(xiàn)滴落流態(tài)。隨著流量比的增大,分散相的慣性力逐漸增大,促進了分散相液線的延長,同時使得分散相的頸部以及斷裂位置向主流道下游方向移動,流態(tài)從滴落態(tài)向擠壓態(tài)和射流態(tài)依次轉(zhuǎn)變。隨著毛細數(shù)的增大,意味著連續(xù)相對分散相的粘性剪切力增大,從而抑制兩相間的界面張力作用,這有助于分散相液線向主流道的下游延伸。受此影響,流態(tài)轉(zhuǎn)變所需的分散相慣性力減小,即流態(tài)轉(zhuǎn)變的臨界流量比隨著毛細數(shù)的增大而減小。

圖10 流態(tài)轉(zhuǎn)變圖Figure 10 Flow regime transition diagram

針對流態(tài)轉(zhuǎn)變圖中的兩條流態(tài)轉(zhuǎn)變線條進行非線性的曲線擬合。得到流態(tài)轉(zhuǎn)變的臨界公式。滴落態(tài)與擠壓態(tài)的臨界關(guān)系如式(11),擠壓態(tài)與射流態(tài)的臨界關(guān)系如式(12):

Ca=0.01Q-1.45,

(11)

Ca=0.13Q-1.99。

(12)

從式(11)和(12)可以發(fā)現(xiàn),擠壓態(tài)向射流態(tài)轉(zhuǎn)變的流量比的相關(guān)參數(shù)的絕對值要大于擠壓態(tài)向滴落態(tài)的,這表明從擠壓態(tài)向射流態(tài)轉(zhuǎn)變所需的流量要大于擠壓態(tài)向滴落態(tài)轉(zhuǎn)變所需的流量。此外,當Ca增大時,流態(tài)轉(zhuǎn)變所需的Q要減小。

4 結(jié) 語

本文在類文丘里十字型結(jié)構(gòu)基礎(chǔ)上,通過改變毛細數(shù)和流量比得到類文丘里十字型流道的3種典型流態(tài)。擠壓態(tài)和滴落態(tài)的液滴形狀分別為柱塞形和圓形。在液滴生成周期內(nèi),兩種流態(tài)的分散相液線長度都隨時間增大,頸部寬度都隨時間先增大后減小。但由于兩種流態(tài)下主導液滴生成的力不同,以及類文丘里結(jié)構(gòu)的影響,使得液滴生成時作用力的競爭關(guān)系發(fā)生改變,造成形貌參數(shù)的變化率呈現(xiàn)出不同規(guī)律。擠壓態(tài)下分散相液線長度變化率表現(xiàn)出波動上升最后趨于平穩(wěn),滴落態(tài)下分散相液線變化率趨勢是先平穩(wěn)不變再迅速上升最后緩慢減小。前者分散相頸部變化率大體表現(xiàn)為先下降后保持不變,后者分散相頸部變化率則是先保持平穩(wěn)再迅速下降最后保持不變。射流態(tài)出現(xiàn)在流量比較大時,由于分散相慣性力足夠大抑制界面張力而無液滴生成。同時,流態(tài)轉(zhuǎn)變的流量比會隨著毛細數(shù)的增大而減小,因此,為保證液滴能夠以滴落態(tài)形式穩(wěn)定生成,應(yīng)適當保證流量比在0～1以內(nèi)。