CFD數值模擬技術在液滴微流控多相流特性研究的應用進展

王炳捷，李輝，楊曉勇，白志山

（華東理工大學機械與動力工程學院，上海，200237）

過程強化是實現綠色工藝的關鍵技術，通過小尺寸和小型化設備的發展，降低單位能耗及副產品生成，并最終達到提高生產效率、降低生產成本、提升安全性和減少環境污染的目的。然而，小尺寸緊湊型設備的發展也給其中流體動力學特性和運動行為的研究帶來了新的挑戰，特別是多相流體系[1-3]。與傳統的宏觀大尺寸系統相比，借助微通道控制的兩相流體具有比表面積大、傳遞距離短、混合速度快等優點，可以減少傳質限制，以獲得更好的性能[4]。因此，液液微流控體系在化學反應[5]、液液萃取[6]、生物分析[7]、結晶過程調控[8]以及結構性材料制備[9]等領域取得了重大進展。為了分析操作條件、流體性質以及通道結構對多相流的影響，當前已經從試驗角度開展了許多有益的研究[2-3,10-13]。

根據以往的試驗研究，微通道內兩相流間的相互作用及動力學行為受到許多參數的影響，如兩相流速、兩相黏度、界面張力、微通道結構及壁面潤濕特性等。同時，預測微通道裝置中液滴尺寸及成型機制的理論分析仍然比較復雜，且準確性不足[14-15]。因此，僅靠試驗研究或理論分析很難對微通道中兩相流動現象和特性有全面的了解。與此同時，微流控系統中參數的測量和流場特性的獲取也是極為困難的工作。由于這兩個原因，目前微流控系統理論多為經驗研究所得。

此外，液滴和氣泡在復雜幾何結構中的運動對于許多科學和工程應用而言非常重要。在過去的幾十年里，無約束介質中氣泡和液滴的運動引起了科學界的極大關注[16]。相反，關于氣泡/液滴在復雜幾何結構通道內演變的研究內容非常有限。一方面微通道壁面可以壓縮氣泡和液滴，在某種程度上，液滴/氣泡的形狀由微通道結構所決定。另一方面，氣泡/液滴可以在無界介質中自由演化，但在微通道結構的限制下運動。為了解決這兩類問題，試驗過程通常將基本現象分離出來，從宏觀上研究第二相的運動。而液滴或氣泡在復雜幾何中演化的問題則通過簡化控制方程的方法來獲得一些近似的解析解[17]。因此，開發能夠實現微通道內液滴成型、液滴聚并/融合、液滴溶解、顆粒聚焦等復雜傳質傳熱過程的數值模擬方法極為必要。

與液滴微流控技術相關的現有文章的綜述重心主要集中在微通道結構和微流控設備的發展[18-19]、微流控技術在獲得結構性/功能性材料的應用[20-21]以及微流控技術及其產物在不同領域的拓展[22-23]，卻很少把目光集中于數值模擬方法對于液滴微流控多相流特性研究的重要輔助作用。本文將從數值模擬所涉及的液滴微流控裝置結構及演變（第2 節）、液滴微流控模擬方法及優化（第3 節）、微通道內兩相流作用過程及原理（第4節）三個方面，系統綜述當前數值模擬方法在液滴微流控技術及應用方面的研究進展（圖1所示），并在第5節針對當前數值模擬方法所存在的瓶頸進行重點闡述，展望未來數值模擬技術在液滴微流控探索應用的發展方向。

1 微通道結構

微通道限制了微流體邊界，因此通道尺寸結構也會對流體特性產生影響。利用黏性剪切力來產生液滴的基礎微通道結構可分為：同軸流（co-flow）、 交 錯 流（cross-flow） 和 聚 焦 流（focusing-flow）。微流控技術的發展也使得微通道結構以上述三種結構為基礎，不斷向復雜組合或優化結構拓展和演變。而針對不同的研究對象（液滴成型、液滴聚并、氣泡溶解等）所采用的微通道結構尺寸和模擬方法以及所考慮的影響因素也不盡相同（見表1和表2）。

1.1 同軸流結構

1.1.1 基礎同軸流結構

用于液-液介質研究的同軸流微通道最早由Umbanhowar 等[53]提出，其結構特征在于分散相流體與連續相流體在平行流中相遇。同軸流結構可以是準二維（2D）平面的或三維（3D）同軸的。其中，前者通常采用多板層組合獲取封閉的同軸流通道結構，而后者則采用不同尺寸的圓形管（或方形管/槽）進行同軸內配合組裝獲得。Lan 等[27]借助CFD 數值模擬方法研究了準二維（2D）平面同軸流結構中的液滴成型過程。通過分析液滴分離位置附近的局部壓力和流速，揭示了滴流成型和射流成型的本質區別［圖2(a)］。Xiong等[26]結合micro-PIV測試和數值模擬發現了液滴形成過程中的內部循環流，并探究了兩相流速、兩相黏度和界面張力對液滴內部循環流特性的影響［圖2(b)］。而Deng等[45]、Guillaument 等[47]和Nacimirad 等[24]則探究了3D 同軸流微流控結構中兩相流速、兩相黏度、界面張力以及微通道壁面潤濕特性等參數對液滴成型方式和尺寸的影響規律。

1.1.2 優化同軸流結構

在傳統同軸流微通道結構的基礎上，通過對分散相或者連續相通道結構優化，可改變微通道內部流體的流動特性。Heuberger 等[31]通過在分散相通道尖端開設方形槽道或布置特定結構的開孔擋板，構建了獨特的噴射式同軸流液滴成型結構，其本身具有雙相流體振蕩器特征，可在高速率和高長徑比條件下連續產生單分散性液滴。隨后的數值模擬幫助確定了上游質量和動量傳遞在激發液-液界面振蕩中的重要作用，而試驗結果則進一步證明了該種優化結構在提升液滴成型的穩定性和產量方面的潛力。而Wang[46]和Chekifi[32]等則改變兩相流體接觸段的幾何結構，構建了錐形同軸流結構微流控系統，并利用數值模擬方法研究了流動參數和錐形結構對液滴成型過程中液滴尺寸和生成頻率的影響。如圖3所示。

圖1 數值模擬方法在液滴微流控技術中的應用

圖2 基礎同軸流微通道結構中液滴成型過程

1.2 交錯流結構

1.2.1 基礎交錯流結構

交錯流結構是指分散相流體和連續相流體以一定角度θ（0°＜θ≤180°）在微通道中相遇[54]。在傳統交錯流結構中，分散相和連續相流體通常是在T形結構通道內以正交方式相遇的。由于T形通道兼具結構簡單性及產物的高單分散性（粒徑偏差CV＜2%）特征，從而被廣泛應用于液滴/氣泡的生成[55]。Yeh等[51]模擬研究了交錯流微通道結構中的多相流流動特性，揭示了單分散海藻酸鈣顆粒形成過程及尺寸調控方法。而Padoin等[43]借助數值模擬方法評估親疏水壁面潤濕特性變化對微通道中氣液流動模式的影響。模擬結果表明，表面張力主導了微通道內部氣液流動特性，具體表現為在相似的邊界條件下，微通道壁面親疏水特性的改變（由親水壁面θ=25°變為疏水壁面θ=105°），造成其中氣液流動模式的改變（由泰勒流變為分層流）。如圖4所示。

表1 模擬方法及對象

1.2.2 優化交錯流結構

微通道的幾何形狀對微通道的動態特性和總壓降有顯著影響，通過改變主通道結構會對其中液滴/氣泡的成型產生顯著的影響。Chandra 等[35]通過在兩相流體交匯處下游主通道內側等間距布置不同高度的三角形阻塞單元，首先構建了優化交錯流結構微流控模型。隨后，采用數值模擬方法分別研究了優化結構對液滴成型、尺寸、形貌和壓降的影響，證明了優化結構對流體傳質傳熱的有益作用。研究結果表明，在所形成的液滴內部及界面處的混合效果顯著提升。此外，三角形阻塞單元間距的變化也會影響液滴成型、液滴長度和壓降。而Zhu等[36]更是通過改變主通道橫截面形狀（梯形、三角形、圓底長方形、半圓形等），構建了優化交錯流結構微流控模型，探究了不同結構對液滴分離時間和尺寸形貌的作用關系。結果表明，對于不同界面的微通道，液滴分離時間、分離直徑和去除時間服從以下規律：三角形截面＜梯形截面＜底部彎曲矩形截面＜矩形截面＜倒立梯形截面。此外，已有報道在微流控裝置周邊布置特殊場（離心場、磁場、電場等），借助外加場作用克服流體原本的物理化學性質（黏性力和表面張力），強化和調控微通道內部流體的流動形式，進而實現流體的輸送、混合等過程。Ren 等[56]模擬研究了離心場作用下并聯交錯流微通道中液滴成型和融合過程，揭示了離心場轉速與液滴尺寸、液滴產率、融合速率之間的構效關系，獲得了液滴成型方式轉變的調控機制。如圖5所示。

表2 微通道結構與特征尺寸

1.2.3 交錯流結構陣列

實際應用中，單通道交錯流微流控結構不能滿足液滴高通量可控合成的要求。考慮到微通道并聯不存在放大效應，因此可以通過并行交錯流結構中的分散相通道，共用連續相通道，構成交錯流結構陣列。在這種結構中，分散相流體被迫同時進入多個交錯流結構微通道，然后在兩相交匯處被連續相流體截斷，同時形成多個擬球型液滴。Kobayashi等[40]采用數值模擬方法探究了交錯流結構陣列中通道尺寸和分散相流速對液滴成型過程的影響，并獲取了該過程中動態界面張力的變化規律。研究表明，最終的液滴尺寸與通道尺寸成正比，而液滴生成速率與通道尺寸成反比。van Dijke等[42]在單分散油滴高通量合成過程的模擬研究中發現，液滴形成點之間的距離由分散相通道末端平臺的高度和施加的壓力決定，而與單分散液滴生成過程中的液滴或氣泡大小無關。如圖6所示。

圖3 優化同軸流微通道結構中液滴成型過程數值模擬［31-32，46］

圖4 基礎交錯流微通道結構中兩相流體相互作用過程數值模擬［38，43，51］

圖5 基礎交錯流微通道結構中兩相流體相互作用過程數值模擬［35-36］

圖6 交錯流結構陣列中兩相流體相互作用過程數值模擬［40，42］

1.3 聚焦流結構

在聚焦流結構微通道中，兩相流體通過流體動力學進行聚焦，然后以延伸流狀態通過收縮的微通道，繼而獲得尺寸更小的液滴。與同軸流結構微流控裝置類似，聚焦流微流控裝置也分為準二維（2D）平面和三維（3D）軸對稱類型。Lan 等[28]在模擬聚焦流結構通道中液滴形成過程時發現，液滴形成過程的總傳質比可達50%。在液滴形成過程的早期，由于連續相的剪切作用，出現了較強的內環流。這增強了兩相之間的傳質，并在液滴頂部形成了一個濃度缺口。在液滴形成過程的后期，液滴運動速度明顯增加，液滴內部不再有循環流動，傳質速度變慢。Dziubiński[52]通過建立準二維（2D）平面聚焦流結構微流控模型，探究了兩側連續相不對等流速條件下，聚焦流流型和位置的變化規律，提出了描述聚焦流形態和位置的簡單理論模型。該模型可在誤差高達25%的情況下確定聚焦流的位置和形狀，有助于推動單顆粒流變學發展，使得借助水動力聚焦來直接控制顆粒受力和運動成為可能。Dang 等[50]在兩相微通道夾角θ為30°的聚焦流微通道結構中，系統研究了通道壁面的接觸角、表面張力和流體黏度對泰勒氣泡長度、體積和形狀的影響。研究發現，在最高表面張力和最高流體黏度的條件下，隨著接觸角的增大，氣泡長度大幅度減小。這是氣泡端部形狀由凸形向凹形的變化和氣泡周圍液膜體積減小的共同作用結果。針對高黏度流體，液滴粒徑控制和成型將更加困難。為解決這一問題，Li等[30]通過對聚焦流結構微流控裝置外加電場作用，模擬研究了電場和電荷在流體界面上的相互作用特性，實現了電場力對高黏度液滴成型動力學的有效調控。如圖7所示。

2 模擬方法

由于微尺度裝置中存在著較大的比表面積，因此，流體界面在許多物理過程中扮演著重要的角色，例如液滴碰撞、合并和破裂[57]。在具有移動邊界和界面流動問題的模擬過程中，存在基于固定網格的離散Navier-Stokes（NS）方程，并以顯式或隱式來表達界面演化的兩種方法，即界面追蹤方法和界面捕捉方法。其中，移動邊界的界面追蹤基于計算網格對界面的顯式描述，具有較高的數值計算精度，但其適用范圍僅限于Stokes流。常用的前沿追蹤法（front-tracking method）[58]屬于界面追蹤方法。與界面追蹤方法不同，界面運動可以簡單地通過界面捕捉方法中對應的相位函數來得到。界面捕獲方法使用單獨的相位函數，以隱式表示界面。用于解決不可壓縮兩相問題的常用界面捕捉方法包括水平集法（level setmethod，LS）[28,59]、流體體積 法 （volume of fluid method， VOF）[60-61]和格子-玻爾茲曼法（lattice Boltzmann method，LB）[62-64]。以上方法對比見表3。

圖7 聚焦流微通道結構中兩相流體相互作用過程數值模擬［28，50，52］

表3 模擬方法對比

2.1 LS方法

2.1.1 保守LS方法

保守LS 方法的關鍵在于附加維度的引入[65]，是基于空間曲面的隱函數表達，最初由Osher 和Sethian[66]于1988 年提出，隨后在1994 年Sussman、Smereka 和Osher 等將其進一步發展。該方法可以簡單地追蹤移動的界面和形狀，精確求解具有復雜拓撲變化的不可壓縮兩相流問題，并準確地表示界面法向和曲率等界面變量[25]。目前已成為模擬不可壓縮多相流問題的一種有效方法，原理如圖8 所示。此外，由于LS 方法采用的是光滑的距離函數來捕捉相界面，各個物理量可以在界面上光滑連續地過渡，相界面的捕捉效果好。Bashir 等[67]采用保守兩相LS 方法模擬交錯流微通道內的液滴成型過程，重點研究了微通道壁面潤濕性對液滴成型的影響。研究結果表明兩相流中液體之間的競爭潤濕特性會導致不穩定的流型，從而破壞均勻液滴的正常成型。特別是對于疏水潤濕區域而言，連續相流體毛細數Cacont小于0.02 時，接觸角的影響作用越來越顯著，后續試驗也驗證了該模擬結果。Naeimirad 等[24]同樣采用保守LS 方法模擬受驅動力影響的不同界面形貌，探究液滴成型方式轉變與系統量綱為1 參數（流速比、黏度比、毛細管數等）變化的構效關系。結果表明，增大流速和降低界面張力會減小液滴尺寸，從而導致液滴成型方式由滴流向噴射轉變。

2.1.2 優化LS方法

采用水平集函數的經典控制方程，在經過多步計算后，流場的畸變將導致水平集函數的梯度變得過大或者過小，造成水平集函數本身的失真。該問題被稱為“質量損失”，是保守LS 方法所固有的，會導致界面的法線、平均曲率和高斯曲率的不準確近似，并且不能通過提高計算精度而解決。為了克服“質量損失”問題，保持水平集函數為一個真正的符號距離函數，需要對它進行重新初始化，此過程被稱作“重定距”[28]。然而，函數初始化過程總伴隨著界面位置的移動，造成質量的損失，導致質量不守恒。此外，改善初始化步驟來矯正質量守恒又會增加計算時間，提升計算成本。近年來，通過優化控制方程、添加新控制方程或耦合其他模型來獲取優化LS方法，已成為解決保守LS方法“質量損失”問題的新途徑。Lan 等[27]通過在經典水平集控制方程中引入兩個額外的附加項，來補償水平集函數的失真，數值模擬所得液滴尺寸的預測值與試驗測量值保持較高的一致性。在隨后的研究中，Xiong 等[26]同樣借鑒該方法來研究同軸流微通道內的液滴成型過程中流場變化和旋流強度。Lan 等[28]通過在水平集函數中設置一個新的控制方程來克服“質量損失”，并將同軸流微通道中液滴形成過程的模擬結果與已有報道的試驗結果相比較，驗證了優化LS 方法的可行性和準確性。而Marchandise等[68-69]則采用不連續的Galerkin（DG）方法來執行水平集輸運方程，確保了質量守恒。Olsson 等[70-71]利用雙曲正切函數代替傳統的符號距離函數，進而減少質量守恒誤差。Wong 等[29]將保守LS 方法與Carreau-Yasuda 應力模型結合，用以監測交錯流通道內液滴成型過程中的尺寸變化。通過經驗獲得的流體流變和物理特性數據服從該應力模型，為模擬計算提供了理論依據。模擬結果表明，剪切成型液滴的尺寸變化由被測流體的流速、黏度等物理特性控制。隨著流速比的增大，數值模擬計算對液滴平均直徑的預測越準確。在測試的最高流速比（0.125）下，水滴平均直徑預測值與測量值的差異約為11%，而在較低的流速比下，則顯示出相當好的一致性。Li 等[30]通過在保守LS 方法中耦合靜電模型，模擬研究了外加靜電場對液滴成型過程影響。數值模擬表明，電場和電荷相互作用產生的電場力，對液滴形成動力學起著重要的控制作用。Gutiérrez等[17]在保守LS方法的基礎上，通過引入拉格朗日-歐拉函數和浸入邊界方法，提出了一種研究復雜幾何形狀中液滴和氣泡演變問題的新方法。具體而言，首先利用保守LS 方法來處理多相域，同時控制質量守恒。隨后引入拉格朗日-歐拉函數來優化仿真區域，在此過程中控制體積的大小被均勻化，網格在重要區域得到細化。因此，移動網格將跟隨氣泡的運動，從而減少了計算域大小、提高了網格質量，顯著降低了計算資源的消耗。最后，利用浸入邊界方法處理復雜的幾何結構，并在拉格朗日-歐拉框架內重新生成內部邊界。該研究方法能夠處理完整的非結構化網格，增加了該模型的適用性。

圖8 LS方法原理

2.2 VOF方法

2.2.1 VOF基礎方法

使用最為廣泛的界面捕捉模型是由Hirt 和Nichols[72]于1981 年首次提出的VOF 方法。該方法采用幾何重建策略來構造流體界面，利用VOF 函數追蹤不同相的體積分數，并采用分段線性界面計算（piecewise-linear interface calculation，PLIC）方法來逼近多相界面[73-74]，原理如圖9 所示。VOF 方法已被用于探究微通道內的液滴成型、聚并、融合、溶解等各種傳熱傳質并存的過程[48]。Padoin等[43]和Kashid等[39]利用VOF方法分別探究了交錯流微通道壁面潤濕特性對等溫氣-液流動模式（泰勒流和分割流）和柱塞狀液滴成型機制的影響規律。Chen 等[37]采用VOF 方法研究了交錯流微通道壁面潤濕特性對液滴成型過程的影響，獲取了毛細管數Ca和通道壁面潤濕角與液滴成型方式和液滴尺寸控制的相互作用關系。Sontti 等[34]采用VOF 方法研究了交錯流微通道中非牛頓流體液滴的成型過程和方式，并獲得了兩相流速、界面張力等因素對液滴成型機制的影響規律。Glatzel等[75]則針對打印機噴頭墨滴生成的實際問題，基于三維模型和VOF 方法模擬研究了墨滴的生成過程及影響因素。

此外，還有眾多學者采用VOF 方法模擬研究特殊結構微通道中的液體動力學行為。針對優化的同軸流微通道結構，Heuberger 等[31]采用VOF 方法研究了分散相通道尖端開設方形槽道的特殊同軸流微通道裝置內部的流體動力學行為以及液滴的成型方式。Chekifi[32]基于VOF方法，研究了在低雷諾數Re和毛細管數Ca條件下，錐形同軸流結構微流控裝置中液滴的生成過程，確定了液滴成型方式、液滴形狀尺寸、產生頻率與流動條件（速度比、黏性效應）和通道幾何尺寸的關系。針對優化的交錯流微通道結構，Chandra 等[35]借助VOF 方法研究了主通道齒形壁面結構（齒形間距、大小等）對通道內部柱栓狀液滴成型、大小、形狀以及齒間壓降的影響，獲取了通道壁面膜和壓降等關鍵特征參數。Zhu等[36]則結合三維模型和VOF 方法模擬研究了主通道橫截面幾何形狀變化對微液滴演化、移動等動力學行為的影響，驗證了液態水對低溫燃料電池中陰極氣體微通道幾何形狀的敏感性。Ren 等[44]利用VOF 方法模擬研究了離心場作用下，并聯交錯流微通道中液滴/氣泡的成型、聚并以及分裂行為。針對交錯流微通道陣列結構，Kobayashi 等[40-41]和van Dijke等[42]均采用VOF方法模擬研究了交錯流結構陣列中多液滴成型過程以及速度場分布和壓力梯度的變化規律。

圖9 VOF方法原理

而在微通道內兩相傳熱傳質方面，Dai 等[33]基于VOF 方法，開發出一個用來分析和預測微通道內部液-液泰勒流的流動特性和傳熱行為的通用模型。同時，借助該模型證實了傳熱系數與流動條件之間的強烈依賴關系，解釋了試驗過程中液-液流動傳熱測量的巨大不確定性和困難。此外，模擬的準確性也被試驗所驗證（大部分Nusselt 試驗測量值都落在數值模擬預測值±10%的范圍以內）。而Li等[38]結合micro-PIV 和VOF 方法，探究了子彈狀液滴的生成過程，預測了液滴的尺寸和流速，獲取了液滴內外流場特性。相關研究發現將有助于優化微通道中柱塞狀/子彈狀液滴結構特性，從而提高多相流傳質傳熱和混合效果。

2.2.2 優化VOF方法

盡管VOF 方法能夠模擬需要追蹤界面的多相問題，但VOF 函數中的體積分數在空間上是一個不連續的階梯函數，造成了兩個相鄰網格的界面也是不連續的，且相關物理量在通過界面時也是不連續的，這個現象稱為寄生流動（parasitic current）。這些非物理速度是由于界面曲率離散化的數值誤差而產生的，常被用來計算表面張力。然而表面張力計算中所帶入的數值誤差被傳遞到動量方程中，此時，則需要速度項來平衡該誤差，進而導致了偽速度。偽速度和寄生流動在微觀尺度上更加嚴重，因為毛細管力占主導地位，這可能會影響微觀尺度上模擬的準確性[76]。目前優化VOF方法的主要途徑就是緩解數值方法造成的偽速度和寄生流動現象[48]。最近，多種用以最小化VOF 模型中的偽速度，進而緩解偽速度問題的創新方法已經被開發，并成功用于研究涉及孔尺度和微流體的微尺度問題[76-79]。

Wang[46]采用VOF 方法和連續表面力（continuum surface force）模型分別捕捉氣-液界面和模擬流體表面張力，對錐形同軸流結構微流控裝置中流體狀態參數和通道外尺寸參數與泰勒氣泡成型方式、頻率和尺寸的構效關系進行了數值研究。類似地，Deng 等[45]同樣采用VOF 模型和連續表面力模型探究了可控尺寸的單分散微液滴的成型過程。模擬結果表明，連續相流體的黏滯阻力、分散相流體的慣性力和液滴界面張力之間的平衡在很大程度上決定了液滴的成型方式及尺寸大小。并在模擬研究的基礎上提出了預測液滴直徑的經驗關系式，同時得到了相關試驗的驗證（液滴直徑的預測值與試驗數據高度吻合，最大誤差為9.2%，平均誤差僅為3.8%）。Guillaument 等[47]使用VOF 模型和單流體（one-fluid）模型組合來追蹤二維微通道中兩相流界面，探究了同軸流微通道結構中界面張力和微通道壁面潤濕特性對超臨界CO2流體成型方式及相位反轉機制的影響規律。Soh 等[48]利用VOF 方法捕捉多相物理量，并在表面張力計算中應用平滑操作以最小化偽速度。同時，耦合多相/多組分模型對不同組分的遷移進行追蹤，模擬研究了交錯流微通道結構中CO2氣泡成型和CO2氣泡在連續相流體硅油中的溶解過程。該優化模型有助于推進微通道中液滴溶解、血管中涉及傳質的藥物運輸以及CO2強化采油率等過程的研究。

2.3 LB方法

2.3.1 LB基礎方法

格子-玻爾茲曼（LB）方法是統計物理學和計算科學的最新發展[49]。該方法是局部的，幾乎不需要進行集群、多線程或GPU 并行處理，與順序處理相比，加速比最高可達20。因此，該方法特別適用于不可壓縮的流動和多物理場（熱流體、電流體動力學、磁流體動力學、傳質和反應流）等流動模擬問題[65]，原理如圖10 所示。與LS 方法不同，LB 方法并不會直接追蹤界面。相反，該模型使用網格來定義整個求解域的節點矩陣，并將流體視為有質量無體積的粒子。粒子在網格的節點間流動，并在節點處相互碰撞。通過模擬確定這些粒子在各個方向上的分布函數，可從這些分布函數中計算得到諸如密度、速度等宏觀參數。該分布函數有4種常用模型：著色模型、Shan-Chen模型、自由能模型和He 模型。著色模型以不同的顏色區分相位，并根據相鄰粒子之間的相互作用對界面進行追蹤。Shan-Chen模型則將多相流體視為非理想流體，可在局部動量不守恒的情況下，很好地解決相變問題。自由能模型使用總密度和密度差來確定兩種流體的密度。He 模型主要用于不可壓縮流體，通過主函數來追蹤界面，其中界面張力被分子間的相互作用所取代。Frisch 等[80]在其研究中證明，可以使用有限差分近似將LB 模型視為Navier-Stokes 方程的穩定形式，對于復雜幾何結構的模擬具有較高的計算效率和精度。Wang等[81]在研究液-固傳質特性的數值模擬過程中發現，使用LB 方法可在更小網格單元尺寸和更短計算時間的條件下，獲得相同精度的計算結果。

圖10 LB方法原理

2.3.2 優化LB方法

然而在眾多模擬研究中發現，宏觀尺度模型在簡單幾何結構中的應用要明顯優于LB模型[82-84]。這是因為，對于LB 方法而言，當使用反彈邊界條件時，模擬的準確性可能會降低。為了獲得較高精度而采用半反彈策略，則會增加計算時間。因此，針對簡單幾何結構，提高效率并保持高精度是LB 方法優化的重要方向[65]。Riaud 等[49]開發了一種模擬兩相微流體相互作用的優化LB 方法。該方法使用Latva-Kokko 再分配策略實現相分離，而界面張力則通過Shan-Chen模型的表達式求解。通過引入粒子碰撞后的再分配步驟，提高算法的穩定性，同時使界面厚度保持恒定，不受界面張力的影響。因此，該優化方法提供了獨立的界面張力控制，受黏度影響很小。與連續表面力模型相比，界面曲率更容易獲取，而偽流動可以保持在與著色模型相當的低水平。隨后，他們采用該模型研究了交錯流微通道結構中氣泡成型過程，并結合以空氣/乙醇為工作介質的T形結構微通道中氣泡成型試驗，通過對界面性質、界面張力以及剪切力與界面張力之間競爭關系的研究，驗證了該模型的有效性。

此外，多相/多組分流體系統在微流體領域中廣泛存在，涉及液滴/氣泡在微通道內的成型、聚并、分裂等過程。由于體系具有較大的比表面積，微通道壁面潤濕特性會對微流體流動行為產生較大的影響。在LB 方法中，相界面問題的處理可以簡單地轉化為對流體粒子間以及流體與微通道壁面間相互作用力的控制[85]。通過耦合多相/多組分模型，則可實現微通道中多組分追蹤，以及流體性質、微通道特征參數等對液滴成型過程作用機制的模擬研究。Dupin 等[62]結合LB 方法（自由能模型）和多相/多組分模型，在低雷諾數Re和低毛細管數Ca條件下，模擬研究了二維交錯流微通道結構中多相流動力學行為，并對流體與微通道壁面相互作用、相界面偽速度等進行調控。Yang等[86]以交錯流微通道結構中液滴成型過程為基礎，將浸入式邊界條件引入LB 方法（Shan-Chen 模型），模擬研究了細胞在封裝過程中的三維旋轉和形變特性。

2.4 LS+VOF耦合方法

LS方法和VOF方法是用以模擬研究復雜界面兩相流相互作用規律及動力學行為最為廣泛的兩種方法。在LS 方法中，通過水平集函數追蹤和捕捉界面。由于水平集函數具有連續性和光滑性特征，可以精確地計算空間梯度。然而，LS方法在保持體積守恒方面存在固有缺陷。相比而言，VOF方法本質上體積守恒，因為該模型中計算和追蹤的是每個單元中特定相的體積分數，而不是界面本身。VOF方法的不足在于空間導數的計算，因為VOF函數本身在界面上不連續。由于不當方法致使表面張力離散化和表面曲率近似化，都會產生偽速度。為了克服LS 方法和VOF 方法的各自缺點，ANSYS FLUENT（ANSYS Inc.，美國）在第14.0 版更新中提供了LS和VOF的耦合方法（CLSVOF）。該方法采用分段線性界面構造（PLIC）的幾何重建來執行重新初始化過程。曲率和界面法向由水平集函數計算獲得，而界面的準確位置則通過平衡每個單元中的體積來調節，從而實現由VOF模型計算獲得體積分數。該方法在重定距水平集函數的同時加強了質量守恒，而流體的表面張力和物理特性則用LS模型相類似的方法計算求解[87]。與VOF 方法相比，CLSVOF 耦合方法可以產生更精確的氣液界面，特別是在氣泡與分散相流體分離階段，可獲得的與試驗更加吻合的結果。Dang等[50]在異構化聚焦流微通道中泰勒氣泡成型過程的數值模擬研究中發現，相較于VOF方法，采用CLSVOF耦合方法模擬所得的前后兩個氣泡的模擬長度差異更小，更加吻合試驗結果）。此外，研究結果還表明，在泰勒氣泡成型過程中，通道壁面接觸角和流體黏度對氣泡的形狀產生了較大影響，表面張力對其的影響則相對較小。

3 研究對象

微通道中微流體間相互作用的結果多樣，借助于CFD數值模擬的方法，通過研究流體間的傳熱傳質過程，可獲取微通道中液滴成型、液滴聚并/融合、液滴溶解以及顆粒聚焦等過程原理和調控機制。

3.1 液滴成型

3.1.1 量綱為1參數

液滴成型是微流控研究的基礎，對于微流體應用至關重要。在液滴成型過程中，從注射泵或者壓力控制器引入的能量部分轉化為界面能，促使液-液界面失穩，進而引發離散的液滴從分散相流體中分離[88]。微流體中可以發生各種性質的流體運動，這通常由相互競爭的物理效應所決定，例如力的平衡。兩相流體相互作用形成液滴過程中液滴本身所受到平衡力可以歸納為兩大類：分離力和附著力。盡管微通道結構復雜、尺寸各異，但液滴成型過程中所涉及的物理學特性和潛在成型機制相似（圖11）。多種力（重力、黏性力、慣性力、界面張力等）相互作用的結果是液滴以不同的流動方式成型，具體而言包括：擠壓成型（squeezing）、滴流成型（dripping）和噴射成型（jetting）[89]。液滴生成過程與界面流動有關，通常涉及3 個主要步驟：首先，分散相和連續相流體在兩相交界處相遇，并形成非混溶界面；隨后，非混溶界面發生大變形，進入不穩定狀態；最后，不穩定界面自發破碎，衰變為相互獨立的液滴。因此，界面張力在液滴成型過程中起著關鍵作用，而與界面張力有關的兩個量綱為1 參數（毛細管數Ca和韋伯數We）在液滴成型過程中尤為重要。

（1）毛細管數Ca毛細管數Ca是黏性力與毛細管壓力的比值，見式(1)。由于慣性力本身不取決于通道尺寸，而通道尺寸的減小將導致黏性力和毛細管壓力的增大以及重力效應的削減。因此，如果通道尺寸足夠小，黏性力和毛細管壓力就會占主導地位。這一原因使得毛細管數Ca成為描述微流體中液滴產生過程最常用的量綱為1參數。在微流體中，毛細管力的取值范圍為10-3～10。

式中，uj為流體流速；μj為流體黏度；γ為兩相間的界面張力。其中下角標j既可代表分散相流體（j=disp），也可代表連續相流體（j=cont）。

（2）韋伯數We盡管對于大多數微流體流動而言，流體的慣性可以忽略。但是對于高流速條件下的射流和非線性氣泡形成以及液滴和氣泡從分散相中的脫離，流體慣性很重要。慣性力和毛細管壓力的相互競爭的結果是產生了韋伯數We。對于大多數微流體流動，We＜1，見式(2)。

式中，ρj為流體密度；d為通道尺寸。

3.1.2 擠壓（squeezing）成型

在毛細管數Ca較小的條件下，此時黏性力被通道壁的限制所取代，液滴以兩個階段的擠壓方式成型，即液滴的出現和生長、與分散相流體分離。在第1階段，分散相流體在連續相流體施加的側向剪切力作用下，在流體流動方向的主通道內不斷積聚。隨著分散相流體前端液滴的不斷增大，液相界面和通道壁面的間隙不斷減小，并在液滴上形成壓力梯度。當此壓力梯度足以克服分散相液滴內部壓力時，分散相界面擠壓變形，形成頸縮。在第2階段，在動態界面張力所驅動的Rayleigh-Plateau 不穩定作用下，分散相流體前端液滴在頸縮處分離，形成長方形或者子彈型（plug-shaped）的柱栓狀液滴/氣泡。擠壓方式產生液滴過程中，液滴在脫落前幾乎完全阻塞兩相交界區域，液滴的大小與毛細管數無關，而由通道的尺寸和兩相流速比決定[34,37,89]。Wang[46]和Sontti 等[34]以擠壓方式在錐形同軸流結構微流控裝置中分別獲得了子彈型的柱栓狀氣泡和液滴。而Chen 等[37]則以擠壓方式在交錯流結構微流控模型中獲得了長方形柱栓狀液滴。

圖11 液滴成型過程模擬

3.1.3 滴流（dripping）成型

當毛細管數Ca繼續增大（更高黏度和/或流速），液滴成型方式由擠壓成型轉變為滴流成型，其中拖曳界面破裂的黏性力占主導地位，而界面張力作用則確保初生的液滴穩定不破裂。在同軸流和聚焦流結構微通道中，液滴滴流成型過程同樣分為兩個階段：液滴的出現和生長；與分散相流體分離。與擠壓方式成型略有不同，由于較大的黏性剪切力，分散相流體前段不斷增大的液滴在完全阻塞微通道之前，即被連續相流體在頸縮處剪斷，因此，所生成的液滴尺寸小于微通道尺寸，保持擬球型形貌，具有高單分散性特征[37]。在這種情況下，液滴尺寸既取決于流速，也取決于兩相黏度比。高黏度的連續相流體在液滴與通道壁面之間將會遇到更大的流體阻力，頸縮位置與分散相流體前沿液滴距離將進一步縮短，其結果是液滴分離過程速度更快，所獲得的液滴尺寸更小。而對于交錯流結構微通道而言，液滴/氣泡在從分散相流體分離后會首先附著在微通道壁面上。因此，存在液滴與微通道壁面分離，并在自身界面張力作用下形成擬球型的自由液滴的第3階段。在此階段，液滴成型的持續時間與通道壁面潤濕特性高度相關。

3.1.4 噴射（jetting）成型

通過增加分散相或者連續相流速，毛細管數Ca將進一步增大，液滴由滴流方式成型向噴射方式成型轉變[90]。在同軸流和聚焦流結構微通道中，液滴噴射成型過程依然分為兩個階段，即液滴的出現和生長、與分散相流體分離。但與擠壓和滴流方式成型均不相同，在液滴噴射成型過程中，由于界面（Rayleigh-Plateau）不穩定性，當連續相流體和分散相流體的慣性力所施加的黏性力大于界面張力時，兩相流體在交匯處形成連續的射流，直觀可見分散相流體具有明顯的“拖尾”現象，此射流在下游末端斷裂成擬球型液滴[37]。由于射流狀態的不穩定性，且液滴從分散相流體脫離過程受兩相黏度和流速等多方面原因的影響，致使噴射方式所獲得液滴具有多分散性特征。此外，液滴形成位置向下游移動，液滴成型位置與液滴之間的距離隨著連續相流體Ca值的增大而增大。當Ca值超過臨界值時，通道內將不能產生液滴，兩相流體流動變為穩定的平行流動。為了避免液滴噴射成型，在高流速下應采用低黏度的連續相。而在交錯流結構微通道中，液滴成型過程仍然包含液滴與微通道壁面分離的第3階段，且分離時間與通道壁面狀態和潤濕特性高度相關。已有眾多學者[27,32,45-46]結合試驗觀測和CFD 模擬方法，系統地研究了液滴/氣泡在不同微通道結構中的滴流和射流成型過程，并從中獲取兩種成型方式的轉換機制和調控方法。而Chen 等[37]和Guillaument等[47]更是探究了微通道壁面潤濕特性等因素對液滴/氣泡成型的影響。

3.1.5 成型方式轉換

隨著毛細管力的增加，分散相流體前端的液滴部分地、暫時地阻塞了連續相流體的流動，從而形成了液滴從擠壓成型到滴流成型的過渡狀態。由于此狀態下，分散相流體局部或間歇地阻塞連續相流體流動，剪切應力和擠壓壓力的共存，因此，過渡狀態下的液滴尺寸既取決于毛細管數Ca，也取決于兩相流速比[37,89]。此外，大量試驗和數值模擬結果表明，3種液滴被動成型方式可以通過改變與界面張力有關的量綱為1 參數（毛細管數Ca和韋伯數We），從而實現相互轉換。然而，相較于分散相毛細管數Cadisp和連續相韋伯數Wecont，連續相的毛細管數Cacont對液滴成型方式有較大的影響，但不同剪切形式所對應的Ca值也存在較大差異（見表4）。這是由于液滴成型方式是由兩相組分、兩相黏度、兩相流速、微通道尺寸結構、邊壁潤濕特性等多方面因素共同作用的結果[91-92]。

3.2 液滴的聚并/融合/溶解

3.2.1 液滴聚并

CFD數值模擬為微通道內液滴之間的碰撞過程提供了更為直觀的研究方法，有助于研究在復雜微通道結構或者多場作用下相同或不同屬性液滴間的相關作用關系。同屬性液滴在微通道中的相互作用過程會發生破碎、聚并和再分裂（圖12）。Anandan 等[93]通過在氣泡邏輯門模型中引入交錯流微通道結構作為氣泡生成通道，模擬研究了氣泡的生成過程，以及氣泡以不同操作條件在兩相鄰微通道切線交匯處的聚并和分裂過程。而Ren等[44]模擬了離心力作用下并聯交錯流微通道中液滴的成型過程，以及液滴在通道交匯處的碰撞、聚并和再分裂過程。結果表明，隨著離心力的增大，液滴成型頻率增大，液滴尺寸減小。對于并聯交錯流微通道而言，微通道與轉動中心間距的增大，致使液滴在成型過程中所受到的離心力作用增大，其結果是距離轉動中心最遠的微通道中所獲得液滴尺寸最小，但液滴成型頻率最快。此外，離心力的作用同樣會加劇并聯微通道所生成的液滴在匯集通道中的碰撞和聚并過程。

3.2.2 液滴融合

不同屬性液滴在微通道中的融合屬于兩相傳質問題，是一個相對復雜的過程（圖13）。Trivedi等[25]借助于試驗觀測和CFD數值模擬研究了交錯流微通道結構中不同屬性液滴間的相互融合過程。研究發現，不同屬性液滴間的融合可以分為接觸、混合、頸縮和脫離4 個連續過程。具體而言，在第1階段，當液滴A被連續相流體輸送到兩相相互作用處時，表面力驅使液滴A界面與分散相流體界面接觸連接，在該過程中毛細作用過程自然發生。在第2階段，液滴A注入分散相中，與分散相流體相互混合，促進了混合液滴B在分散相流體前沿的出現和長大。在第3階段，受到連續相流體施加的橫向剪切力，分散相流體不斷向前端積聚，致使分散相前端液滴不斷長大，分散相后端出現“頸縮”。在第4階段，當液滴增長到臨界尺寸，混合液滴B從分散相頸縮處斷裂，并逐漸脫離微通道頂端壁面，形成自由液滴B。CFD模擬研究表明，在適當的操作條件下，該過程可以自動同步。

表4 液滴/氣泡成型方式及轉換條件

3.2.3 氣泡溶解

微通道中氣泡在連續相流體中的溶解過程同樣屬于兩相傳質問題，試驗方法對于該過程的研究往往存在極大的難度，僅僅停留在針對氣泡尺寸變化的直觀測量。而借助CFD 模擬方法，通過研究氣泡在連續相流體中溶解階段的內外流場變化特性，從而獲取氣泡微通道溶解過程機理和調控方法（圖14）。Soh 等[48]通過耦合傳質模型和VOF 模型，在二維微通道中模擬研究了微通道內CO2氣體剪切成型和硅油中的溶解過程，獲取了CO2氣泡溶解過程中長度和體積的變化規律，導出了溶解氣泡長度和體積隨時間變化的解析解，并與恒定傳質系數模型等3 種模型的模擬結果以及試驗結果進行比較，進而證明了模型的準確性和可靠性。在此過程中，流體體積模型用于捕獲多相物理量，并在表面張力計算中應用平滑操作以最小化偽速度。對于相遷移的多組分追蹤，采用了α因子和驅逐操作方法，以確保在相應的階段追蹤正確的對象。這項研究所開發的耦合模型可用于研究微通道中液滴的溶解、血管內藥物的運輸等傳質過程。

3.3 顆粒聚焦

水動力聚焦是微流體中應用最廣泛的技術之一，已被廣泛應用于化學/生物分析，包括流式細胞術、單分子檢測以及用于快速化學和酶動力學研究的層流混合器[94]。通過在彎曲微流體通道中引入Dean流，Mao等[95]基于三維流體動力聚焦技術，開發了“微流控漂移”技術，實現了單層平板微流體裝置中細胞群或顆粒群的三維流體動力聚焦（圖15）。借助于CFD模擬，顆粒群在微通道內的三維流動聚焦過程得以研究，具體分為3 個階段：第1階段，包含顆粒的流體（以下簡稱為“顆粒流”）和鞘層流分別從平行進口引入微流控芯片中，并在進入90°曲線通道前進行混合；在第2 階段，在曲線通道處，所引入的Dean 流在垂直方向上以雙渦的形式對顆粒進行橫向加速，將顆粒從通道頂部和底部掃向通道中心平面，并進一步將顆粒整體拖向通道中心區域，這一步被稱為“微流控漂移”，有效地實現了顆粒在垂直方向上的流動聚焦；在第3階段，顆粒流受到垂直于流動方向的兩側對稱旁路鞘層流的水平剪切作用，促使原本聚焦在微通道中心區域的顆粒進一步在中心軸線聚焦，實現單顆粒在中心軸線方向的依次排列。該研究成果有助于實現細胞在微通道中的流動聚焦檢測。

圖12 液滴聚并過程模擬［44，93］

圖13 液滴融合過程模擬［25］

圖14 氣泡溶解過程模擬［48］

4 結語與展望

在受限微流控通道中操作兩相流體獲得液滴，進而實現微通道液滴聚并/融合、液滴溶解、顆粒聚焦等復雜傳質傳熱過程，已被廣泛應用于許多科學領域。當下數值模擬方法的快速發展為預測和分析微流控通道內的多相流問題提供了更為直觀有效的幫助，但是隨著微通道結構的復雜化、多相流體系的復雜化以及多場耦合微通道內多相流相互作用的復雜化，也給數值模擬方法帶來了前所未有的挑戰，同時促使著更為完備的數值模擬技術的革新。在未來，數值模擬技術在液滴微流控探索應用的發展方向將會集中于以下幾個方面。

圖15 顆粒聚焦過程模擬［94-95］

（1）數值模擬方法的準確性和適用性的提升。隨著微通道結構復雜化和小型化集成，現有基礎數值模擬方法，包括LS 方法、VOF 方法和LB 方法等，由于本身函數/模型的局限性或缺陷，難以避免試驗結果與數值模擬結果的吻合性較差。此前，一些研究人員試圖通過耦合其他模型或優化本構方程的方式，克服LS方法的“質量損失”、VOF方法的偽速度和寄生流動以及LB 方法的（半）反彈邊界問題，進而提升數值模擬結果的準確性，但涉及組合方法或者優化方法的適用范圍仍需進一步擴展。

（2）模擬對象從牛頓流體體系向非牛頓流體體系的拓展。微通道中液滴成型的研究大多采用傳統的牛頓流體體系，而對非牛頓流體的研究很少。非牛頓流體具有剪切力和變形速率非線性相關的特性，不同的力，包括表面張力和黏性力，控制著液滴動力學。當涉及非牛頓流體時，流體的流變學可能很重要。然而，當前沒有建立單獨的本構方程來描述所有流體的流變圖。了解微流控通道中的非牛頓流體液滴成型的動力學機制對于確保液滴尺寸、形態和產量等特性要求至關重要。此前，一些研究人員專注于分析具有黏彈性特性的液滴的動力學和形變，但對于其他非牛頓流體體系產生可控液滴背后的相關物理學性質的深入了解仍然有限。

（3）液滴輸運過程動力學行為及內部流場特性研究的深入。微通道受限空間內的多相流相互作用過程以及液滴/氣泡在復雜幾何結構的動力學行為對于許多科學和工程應用而言非常重要。此前，這兩方面的研究常常通過簡化控制方程的方式來尋求一些近似的解析解。隨著光學觀測技術日新月異的進步，為揭示復雜化/小型化微通道內多相流相互作用過程中，液滴內部流場特性及變化規律提供了有利的支撐。而數值模擬軟件的發展為定量化預測液滴在受限空間內的動力學行為（自轉、公轉、跳躍）提供了可能，加深了對微通道內傳質傳熱過程強化的認知。