夏侯彬,蔣志韜,陳忠利,錢麗娟,2
(1.中國計量大學(xué) 機(jī)電工程學(xué)院,浙江 杭州 310018;2.浙江省智能制造質(zhì)量大數(shù)據(jù)溯源與應(yīng)用重點實驗室,浙江 杭州 310018)
微液滴以其尺寸小,可控性強(qiáng),比表面積大等優(yōu)點被廣泛應(yīng)用于生物醫(yī)學(xué)[1]、藥物篩選和輸送[2]、食品工業(yè)[3]和化妝品制造[4]等行業(yè)。通常情況下,利用微流控技術(shù)生成微液滴主要依賴于多相流的流體流動,但同時也與所選擇的微流道幾何形狀密切相關(guān)。這是由于改變微流道的幾何形狀可以有效地控制其內(nèi)部多相流的流體動力學(xué)和界面演化過程,從而更好地制備所需特性的微液滴。
科研人員采用實驗和數(shù)值模擬方法,就十字型流道結(jié)構(gòu)變化對液滴生成的影響做了一系列研究。Liu和Zhang利用格子玻爾茲曼方法(lattice Boltzmann method, LBM)對十字型微流道進(jìn)行了三維數(shù)值模擬,發(fā)現(xiàn)當(dāng)微流道的寬深比以及兩相進(jìn)口的寬度比越大時,生成的液滴尺寸也就越大[5]。Gulati等通過實驗探究得出十字交叉處的圓形結(jié)構(gòu)會對液滴生成產(chǎn)生阻礙作用[6]。Zhang等利用水平集(level set method, LS)方法,比較了十字交叉型出口處4種不同頸縮結(jié)構(gòu)下液滴生成的尺寸和周期,發(fā)現(xiàn)當(dāng)出口頸縮結(jié)構(gòu)為梯形時生成的液滴尺寸最小,頻率最高[7]。Ngo和Iqbal等采用流體體積法(volume of fluid method,VOF)對十字流道進(jìn)口夾角進(jìn)行了研究,觀測得出液滴尺寸在夾角為90°時最大,而隨著進(jìn)口夾角偏離90°,液滴尺寸會逐漸減小[8-9]。
綜上所述,十字型流道結(jié)構(gòu)變化的研究多集中在結(jié)構(gòu)變化對于液滴尺寸和生成周期的影響。然而,對于流道結(jié)構(gòu)變化下液滴流態(tài)的分析以及液滴生成時壓力積聚現(xiàn)象并未有太多提及。查閱相關(guān)文獻(xiàn)發(fā)現(xiàn),文丘里結(jié)構(gòu)可以有效緩解流道內(nèi)的壓力積聚問題,已有前人將文丘里結(jié)構(gòu)加入T型流道中,來減緩液滴生成時的壓力積聚問題[10]。因此,本文采用VOF方法,對類文丘里十字型流道的液滴生成流態(tài)進(jìn)行分析,以便更好地理解液滴生成的機(jī)制,從而為后續(xù)實現(xiàn)液滴精確控制提供一定指導(dǎo)。
VOF界面追蹤方法[11]利用體積函數(shù)α表示某一項流體在一個單元網(wǎng)格中所占體積分?jǐn)?shù)。本文以分散相體積為主要研究對象,每個網(wǎng)格單元中兩相流體的體積分?jǐn)?shù)關(guān)系如下:
αd+αc=1。
(1)
式(1)中下標(biāo)d和c分別代表分散相流體和連續(xù)相流體。
通過求解體積函數(shù)的輸運(yùn)方程來追蹤兩相界面,如式(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ù),可以通過兩相流體的體積分?jǐn)?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)
當(dāng)相界面與通道內(nèi)壁接觸時,需要考慮壁面的粘附效應(yīng)[13],它通常由壁面接觸角決定,且接觸角用于計算壁面附近的界面曲率,如式(10):
(10)
本文所用的十字型微流道的二維數(shù)值模擬模型如圖1,兩相入口寬度和長度均相等,通道寬度設(shè)置為Wc=Wd=100 μm,通道長度設(shè)置為Lc=Ld=100 μm。將交匯處改為類文丘里結(jié)構(gòu),其中縮進(jìn)長度設(shè)置為Wv=Hv=25 μm,縮進(jìn)角度設(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
為了簡化計算模型,對計算模型進(jìn)行如下假設(shè):兩相流動均為不可壓縮流動;與外界無熱量交換;分散相與連續(xù)相在流動過程中無化學(xué)反應(yīng)發(fā)生;兩相流體的物性參數(shù)為常數(shù)。由于微流道流動中的雷諾數(shù)小于0.1,故選擇層流模型和基于壓力基的非穩(wěn)態(tài)求解器。邊界條件設(shè)置如下:兩相入口都設(shè)置為恒定的速度入口,出口設(shè)置為壓力出口,出口壓力Pout=0 Pa。
使用VOF框架中的有限體積法求解控制方程以及壓力隱式分裂算子(PISO)求解連續(xù)性方程和動量方程。壓力方程采用PRESTIO!算法離散,而動量方程則通過二階迎風(fēng)格式離散。相界面重構(gòu)采用分段線性界面法[14]。在全局庫朗數(shù)小于0.2的標(biāo)準(zhǔn)下,采用變時間步長法確保計算的收斂和穩(wěn)定。松弛因子設(shè)定如下:壓力0.2,密度0.3,體積力0.3,動量0.2。
為了尋找合適的網(wǎng)格尺寸,本章采用結(jié)構(gòu)化四邊形網(wǎng)格對微流道的二維幾何模型進(jìn)行劃分,網(wǎng)格劃分如圖2。選取三種不同網(wǎng)格數(shù)進(jìn)行網(wǎng)格獨立性驗證,即32 400、41 560和64 750。定義兩個形貌參數(shù)來進(jìn)行獨立性驗證,分別是分散相的液線長度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é)果逐漸收斂。當(dāng)網(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
為了驗證模型的正確性,選取Soroor[15]所做實驗進(jìn)行對應(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ù)值模擬中某些時刻液滴的輪廓存在差異,但整體曲率輪廓和文獻(xiàn)中實驗所得的曲率輪廓幾乎吻合,故可以保證數(shù)值模型的正確性。
圖4 數(shù)值模擬和實驗所得的液滴生成過程的形貌演化對比圖Figure 4 Morphological evolution of droplet generation process obtained from numerical simulation and experiment
數(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ù)相的擠壓作用增強(qiáng),分散相的徑向?qū)挾乳_始減小。分散相頭部進(jìn)入主流道,形成柱塞狀。從流場變化可知,連續(xù)相的流線從交匯處的出口轉(zhuǎn)向進(jìn)口,流線的聚集程度減緩,使得交匯處的壓力減小。
頸縮階段(見圖5,t=6.5~10.0 ms):分散相的頸部內(nèi)凹變形并逐漸拉長變細(xì),以抵抗連續(xù)相的擠壓作用。連續(xù)相的流線垂直于分散性的流線,促進(jìn)頸部的凹陷,使得頸部處的壓力增大。
斷裂階段(見圖5,t=10.0~11.5 ms):連續(xù)相流線的垂直分布表面連續(xù)相對頸部的擠壓作用仍然存在,使得壓力繼續(xù)增大,促進(jìn)頸部快速斷裂形成液滴。在表明張力作用下,液滴尾部與斷裂后的分散相前端后退回縮。
為了更好的理解液滴生成的機(jī)理,采用1.4節(jié)中定義的分散相的液線長度l和頸部寬度δ以及它們各自的變化率(Δl/Δt,Δδ/Δt)定量描述液滴生成過程中相界面的形貌演化。此外,將分散相頭部前端界面內(nèi)外的拉普拉斯壓差(Δp=pi-po),見圖5(a),t=0 ms,隨時間的變化作為考察對象。從圖5可以看出,在膨脹階段,分散相頭部在兩相交匯處體積持續(xù)擴(kuò)大,l和δ逐漸增大,見圖6(a)(c),t=0~3.1 ms。該階段內(nèi)來自側(cè)通道的連續(xù)相對交匯處的分散相的擠壓作用增強(qiáng),進(jìn)而抑制了分散相頭部在y方向上的膨脹,因此Δδ/Δt減小,見圖6(d),t=0~3.1 ms。而這又促進(jìn)了分散相頭部向主流道下游的移動,使得Δl/Δt增大,見圖6(b),t=0~3.1 ms。而后,當(dāng)分散相頭部到達(dá)兩相交匯處出口后,由于流出空間的減小,分散相在y方向?qū)B續(xù)相的擠壓作用的抵抗增強(qiáng),故Δ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)定性增加,進(jìn)而液滴前端界面趨向圓形。而在膨脹階段后期,連續(xù)相的擠壓作用加上流出空間的減小,使得分散相頭部前端界面不斷前凸,界面曲率半徑減小,進(jìn)而Δp增大,見圖6(e),t=0~3.1 ms,液滴穩(wěn)定性減弱,液滴前端界面偏離圓形。
在擠壓階段,隨著連續(xù)相擠壓作用的持續(xù)增強(qiáng),交匯處內(nèi)凸起的分散相受到擠壓,δ開始縮小,Δδ/Δt變?yōu)樨?fù)值且其絕對值隨著時間持續(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)定性進(jìn)一步減弱,液滴前端界面由圓形變?yōu)榧忓N形。
從擠壓階段過渡到頸縮階段時,兩相交匯處被分散相完全堵塞,連續(xù)相的擠壓作用達(dá)到最大,頸部開始內(nèi)凹,Δδ/Δt達(dá)到最小,見圖6(c),t=6.5~10.0 ms,并趨于穩(wěn)定。而且,隨著頸部收縮程度和分散相頭部體積的增大,頸縮過程對l的促進(jìn)作用大幅降低。因此,該階段內(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]的誘導(dǎo)下斷裂,δ減小到0,由于分散相對出流口的堵塞,Δδ/Δt并無較大變化,見圖6(c)和(d)。由于斷裂處頸部的體積很小,液滴頸部的夾斷對分散相液線向下游移動的影響極其有限,因此,在該階段內(nèi),分散相的液線長度延續(xù)了頸縮階段結(jié)束時的變化率,即Δl/Δt和Δδ/Δt保持不變,見圖6(b)。此外分散相頭部前端的界面形貌已經(jīng)達(dá)到穩(wěn)定,故Δp基本不變,見圖6(e),液滴前端界面保持圓形不變。
在擠壓態(tài)的基礎(chǔ)上,減小流量比和毛細(xì)數(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)液滴流體動力學(xué)行為的演化過程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)向進(jìn)口,聚集程度減緩,交匯處的壓力略微減小。
頸縮階段(見圖7,t=18.50~21.50 ms):頸部受到連續(xù)相粘性剪切作用開始內(nèi)凹變薄,主流道內(nèi)分散相體積緩慢增大。連續(xù)相流線垂直于頸部,流線匯聚使得壓強(qiáng)升高。
斷裂階段(見圖7,t=21.50~23.75 ms):在連續(xù)相粘性剪切作用下,流線匯聚位置相較于擠壓態(tài)更靠近分散相入口,進(jìn)而交匯處壓強(qiáng)繼續(xù)增大,促進(jìn)頸部斷裂形成液滴。
滴落態(tài)液滴界面演化特征參數(shù)如圖8。在膨脹階段,連續(xù)相的粘性剪切作用作為主導(dǎo)[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
進(jìn)入擠壓階段,連續(xù)相的粘性剪切作用持續(xù),l和δ繼續(xù)保持穩(wěn)定增長,Δl/Δt和Δδ/Δt保持不變,見圖8(a)~(d),t=13.25~18.50 ms。在擠壓階段后期,由于分散相頭部體積的增大,使得連續(xù)相的粘性剪切與拖曳作用持續(xù)增強(qiáng),進(jìn)而引起頸部變細(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ù)相的粘性剪切和拖曳作用的增強(qiáng)速率也隨之降低并趨于穩(wěn)定,進(jìn)而引起Δl/Δt和Δδ/Δt的絕對值的增大逐漸變緩并趨于穩(wěn)定,見圖8(b)和(d),t=18.50~21.50 ms。分散相頭部尺寸的增大,使得前端的曲率半徑繼續(xù)增大,進(jìn)而Δp減小,見圖8(e),t=18.50~21.50 ms。
斷裂階段,Rayleigh-Plateau不穩(wěn)定性引起液滴頸部斷裂,δ變?yōu)?,Δδ/Δt的絕對值急劇增大,見圖8(c)和(d)。不同于擠壓態(tài),盡管液滴頭部體積達(dá)到穩(wěn)定,但此時分散相頭部前端在界面張力的作用下逐漸變?yōu)榍蛐?曲率半徑不斷增大,Δp繼續(xù)減小,見圖8(e)。
在擠壓態(tài)的基礎(chǔ)上,保持毛細(xì)數(shù)不變繼續(xù)增大流量比,液滴流態(tài)轉(zhuǎn)變?yōu)樯淞鲬B(tài)。相較于擠壓態(tài)和滴落態(tài),射流態(tài)下并無液滴生成,而是形成一條均勻穩(wěn)定的液線,隨著兩側(cè)通道的連續(xù)相一起沿主流道流出,如圖9。在射流態(tài)下,分散相的慣性力主導(dǎo)分散相液線的流動,且足夠大的分散相慣性力可以抑制在兩相交界處由界面張力誘發(fā)的Rayleigh-Plateau不穩(wěn)定性,因此分散相與連續(xù)相界面處不會發(fā)生波動和破裂。
圖9 射流態(tài)流體動力學(xué)行為的演化過程Figure 9 Evolution of the fluid dynamics behavior of threading regime
如上所述,毛細(xì)數(shù)和流量比的變化改變了液滴生成時不同力之間的競爭關(guān)系。因此,基于毛細(xì)數(shù)和流量比,繪制了流態(tài)分布圖,如圖10。當(dāng)毛細(xì)數(shù)和流量比都較小時,連續(xù)相的粘性力和分散相的慣性力也很小,分散相在界面張力的主導(dǎo)下在兩相交匯處經(jīng)歷膨脹和擠壓過程,即出現(xiàn)滴落流態(tài)。隨著流量比的增大,分散相的慣性力逐漸增大,促進(jìn)了分散相液線的延長,同時使得分散相的頸部以及斷裂位置向主流道下游方向移動,流態(tài)從滴落態(tài)向擠壓態(tài)和射流態(tài)依次轉(zhuǎn)變。隨著毛細(xì)數(shù)的增大,意味著連續(xù)相對分散相的粘性剪切力增大,從而抑制兩相間的界面張力作用,這有助于分散相液線向主流道的下游延伸。受此影響,流態(tài)轉(zhuǎn)變所需的分散相慣性力減小,即流態(tài)轉(zhuǎn)變的臨界流量比隨著毛細(xì)數(shù)的增大而減小。
圖10 流態(tài)轉(zhuǎn)變圖Figure 10 Flow regime transition diagram
針對流態(tài)轉(zhuǎn)變圖中的兩條流態(tài)轉(zhuǎn)變線條進(jì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)變所需的流量。此外,當(dāng)Ca增大時,流態(tài)轉(zhuǎn)變所需的Q要減小。
本文在類文丘里十字型結(jié)構(gòu)基礎(chǔ)上,通過改變毛細(xì)數(shù)和流量比得到類文丘里十字型流道的3種典型流態(tài)。擠壓態(tài)和滴落態(tài)的液滴形狀分別為柱塞形和圓形。在液滴生成周期內(nèi),兩種流態(tài)的分散相液線長度都隨時間增大,頸部寬度都隨時間先增大后減小。但由于兩種流態(tài)下主導(dǎo)液滴生成的力不同,以及類文丘里結(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)變的流量比會隨著毛細(xì)數(shù)的增大而減小,因此,為保證液滴能夠以滴落態(tài)形式穩(wěn)定生成,應(yīng)適當(dāng)保證流量比在0~1以內(nèi)。