不混溶液體表面上蒸發液滴的動力學特性

李春曦，莊立宇，施智賢，葉學民

（華北電力大學動力工程系，河北保定071003）

引 言

近年來，一種液滴在另一種不混溶液體表面上的鋪展受到了廣泛關注，并在各個領域得到了應用，如石油回收、支氣管藥物的霧化輸送、灌注液體型超滑表面(slippery liquid-infused porous surfaces，SLIPS)的相變換熱等[1-5]。了解蒸發液滴在不混溶液體表面上的動力學特性對上述應用具有重要意義。

鋪展系數通常用于分析液滴的潤濕行為和幾何結構。Harkins[3]提出了一種液體在另一種不混溶液體基底上的鋪展系數，定義如下:

式中，σLi為純液體i的表面張力；σLij為液體i和j 之間的界面張力，1、2 分別表示液體基底和蒸發液滴。當液滴放置于另一種液體上時，是否會自發鋪展形成單分子薄膜還是會去潤濕形成穩定的液滴取決于鋪展系數S 的正負。正鋪展系數(S ＞0)表明有足夠的力驅動液滴自發鋪展并潤濕液體基底[4]，而負鋪展系數(S ＜0)則會導致一種液體在另一種液體表面形成穩定形狀的液滴。圖1(a)顯示了部分潤濕(partial wetting)的情形下，上層液滴的接觸角不為零且穩定漂浮在下層液體上，稱為液體透鏡；圖1(b)為假部分潤濕(pseudo-partial wetting)狀態下，上層液滴形成了一個接觸角為有限值的球形帽(有學者將這種狀態也稱為液體透鏡[5])，其邊緣在底層液體表面鋪展形成了厚度為分子水平的薄膜；圖1(c)顯示的是在完全潤濕(total wetting)的系統中，上層液滴表面平坦且完全覆蓋下層液體，類似于一個扁平的“煎餅”。Sebilleau[5]指出還存在一種準完全潤濕狀態(pseudo-total wetting)，其圖形與圖1(a)近似，不同之處在于液體2不僅形成穩定形狀的液滴漂浮于液體1之上，而且其表面完全被液體1覆蓋。

蒸發對于不混溶液體表面上液滴的運動特性有重要影響，這方面研究以實驗居多。Nosoko 等[6]研究了施加溫度梯度的去離子水表面上揮發性正戊烷透鏡的蒸發特性，發現透鏡純度和溫差對蒸發時間和透鏡幾何結構有顯著影響。Rahman等[7]研究了黏性液滴在液-液界面上鋪展早期的動力學，發現由于液滴底部出現變形，導致鋪展半徑與時間的關系不再符合指數律。Shabani等[8]觀察到油氣界面漂浮水滴的兩種穩定結構：非凝聚液滴和帽/珠液滴。Sun等[9]、Liu等[10]對不同揮發性物質的液滴在水表面的鋪展行為進行了實驗研究，他們均發現液滴的形狀變化經歷了鋪展主導階段(半徑增大)和蒸發(半徑收縮)主導階段。Dussaud等[11]研究了揮發性油透鏡在密度較大液體基底上自發鋪展的動力學特性，指出揮發性的、不混溶液滴的鋪展半徑R與時間t 的關系為R - t1/2，其鋪展指數1/2 小于非揮發性液滴的3/4。Buffone 等[12]對蒸發或非蒸發液體基底上蒸發液體透鏡的形成、穩定性和溫度分布進行了實驗研究，揭示了透鏡與液體基底間的相互作用。Shimizu 等[13]研究了壓力對正戊烷和R113 液滴在水中蒸發的影響，認為蒸發液滴由氣相和未蒸發的液相組成，指出正戊烷和R113氣泡在上升過程中會出現形狀振蕩和軌跡擺動。

近年來，也有很多學者對不混溶液體表面上液滴的鋪展過程展開了理論研究，但現有的理論模型較少涉及蒸發。對液體透鏡的其他研究主要集中在鋪展/去潤濕動力學特性[14-16]、線張力[17]、潤濕轉變[18-19]和可變形基底上的潤濕行為等[20]。Burton等[21]根據Young-Laplace 方程開發了一個分析透鏡幾何形態的代碼，通過求解界面接觸線上的Young-Laplace 微分方程，計算出三個界面的軸對稱輪廓。Craster 等[1]運用潤滑理論建立了透鏡在各種不同參數組合下(包括兩種液體的密度比、黏度比和表面張力比)的鋪展，計算了拉伸應力在透鏡中的作用，同時也為后續的研究奠定了基礎，但其所建模型未涉及蒸發。Greco 等[22]研究了懸浮在兩種流體界面上的液滴在平行于界面的溫度梯度作用下的熱毛細驅動運動，指出界面液滴的遷移方向與溫度梯度相反。George 等[23]研究了在不滿足Neumann 三角的系統中較重液滴在較輕液體上的漂浮，并將這種現象歸因于線張力的作用，它可防止液滴完全浸沒于液體基底。Phan等[24]證明了較重的液體透鏡可以漂浮在較輕的液體基底上。Guan 等[25]發現在SLIPS 表面上，小固著液滴(直徑約為2 mm)的蒸發遵循理想的恒定接觸角模式。Brochard-Wyart 等[26]研究了黏性液滴在非黏性液體上的鋪展情況，并確定液滴鋪展遵循R - t1/4。Gelderblom 等[27]研究了液滴在液體基底上蒸發的情況，對兩種液體的速度場進行了參數化研究，揭示了兩種液體之間的接觸角及黏度比對速度場的影響。另外，在微觀尺度，采用分子動力學模擬的方法對此類問題進行的研究成果也有報道[28]。

圖1 三流體系統可能的潤濕情況Fig.1 Three possible wetting conditions for a three fluid system

綜上所述，目前關于液滴在不混溶液體表面鋪展的研究大多偏重于實驗方向，但其結果受實驗條件所限，存在許多不一致的地方，也不便于深入了解相關機理，而現有的理論模型較少涉及蒸發。為了更好地理解液滴在液體表面的蒸發機理，以促進相關技術的開發和利用，有必要采用理論建模和數值計算的方法對此類現象進行全面的研究。本文采用潤滑理論建立蒸發液滴鋪展的理論模型，通過數值計算，以接觸線的運動等為著眼點，對不混溶液體表面上蒸發液滴的動力學特性進行探究，以增進對其物理機制的了解。

1 理論模型

1.1 物理模型

圖2 液滴放置于均勻加熱的液體基底上的示意圖Fig.2 Sketch of a drop resting over uniform heated liquid substrate

液體i(i=1,2)鋪展和蒸發過程滿足連續方程、動量方程和能量方程：

式中，T*w為固體壁面溫度。

在液-液界面處，滿足速度場、溫度場連續性，切向和法向應力平衡，運動邊界條件以及熱流量連續性：

潤濕性液滴會鋪展成接近分子厚度的薄膜，而非潤濕性液滴則會形成穩定的透鏡形態。本文主要研究非潤濕揮發性的純液滴在不混溶液體表面的鋪展蒸發過程。當界面溫度不均勻時，液膜的表面張力發生變化，產生Marangoni效應并導致液膜產生流動和明顯變形。上層液滴位于h1≤z ≤h2和-∞≤x ≤+∞約束的區域，為了避免三相接觸線處的應力集中，仿照Craster 等[1]的做法，假定上層液滴可在下層液體表面鋪展成一層極薄的前驅膜，即圖1(b)的情形，此時式(1)所描述的三種界面張力模型可簡化為兩種界面張力模型，即液1 與液2 以及液2與空氣的界面張力。該前驅層厚度h∞由分離壓維持，分離壓Π表示為[1]：

式中，h ≡h2- h1；h∞為前驅膜厚度；n ＞m ＞1，指數(n,m)取(4,3)。式(26)中右側第一項和第二項分別代表分子間斥力和引力對分離壓的不同貢獻。

1.2 初始條件和邊界條件

采用PDECOL程序模擬上述演化方程。由于演化方程中的變量及物性參數均為無量綱量，因此除非特別說明，下文給出的表達式及變量分布曲線圖均為無量綱形式。選取-5 ≤x ≤5 的一維區域進行計算，將該區域劃分為均勻網格，并進行網格無關性驗證，在滿足精度和減少計算時間的條件下設定網格數為2000個。

計算初始條件為[1]：

圖3 初始時刻界面高度分布Fig.3 Initial conditions of the interface height distribution

計算中PDECOL程序的邊界條件設置如下：

在加熱條件下，以演化方程式(23)和式(24)為基礎，計算存在蒸發情形時的上層液滴在下層液體上的演化過程。結合相關資料[9,11]可得有量綱基本參數的取值范圍(表1)，并由此推算出模擬中所用無量綱參數的取值(表2)。模擬中涉及的參數取值如下：Ca1= Ca2= 0.1，Bo = 10-1，0.2 ≤σ ≤5，0.25 ≤μ ≤10，0.1 ≤ρ ≤10，γ = 5，λ= 1，E = 0.001，Σi= 0.001，h1∞= 1，h∞= 0.01。

表1 典型參數的取值或變化范圍Table 1 Order of magnitude estimates for the relevant physical parameters

表2 無量綱參數定義及其數量級Table 2 Representation and order of magnitude estimates of the relevant dimensionless groups

2 結果與討論

2.1 演化過程

圖4 顯示了不混溶液體表面上液滴蒸發的時間-空間演化過程，圖中液滴下半部分曲線對應表示液-液界面h1，上半部分曲線對應液-氣界面h2，h2和h1之間的部分即為由液體2 形成的液滴（其邊緣存在一層厚度極薄的預置液膜），后面的圖也是與此相同的表示方法。液滴在蒸發過程中接觸線后移，液滴厚度和體積逐漸減小。由式(19)和式(20)可知在加熱條件下，液滴界面溫度與液滴厚度呈負相關。蒸發過程中液滴下表面與液滴上表面溫度差引起表面張力梯度差，產生的Marangoni 效應使液體分子從低表面張力處流向高表面張力處，即從熱區域流向冷區域。由于蒸發對液滴變形的影響，從液-氣界面蒸發去除的液體分子必須得到補充以滿足界面力平衡，導致液滴內部產生從下向上的質量遷移，并達到瞬時準平衡態[9]。隨著蒸發的進行，更多的液體分子進入周圍的空氣中，因此液滴上表面的液膜高度隨時間降低，下表面的液膜高隨時間先降低后升高，液滴體積不斷減小，最終完全蒸發。蒸發、重力、毛細力和Marangoni 效應的耦合作用是導致蒸發液滴出現此種現象的主要因素。

圖4 液滴快速蒸發到最終形狀的時間-空間演化Fig.4 Spatial-temporal evolution of the fast evaporation of a drop to a final shape

從圖5可以看出，演化過程分為兩個階段：①t =0～0.1，“鋪展主導”階段，在該階段內，液滴在毛細力和重力作用下緩沖到液體基底中，初始條件的影響被消除；②t = 0.1～100，“蒸發主導”階段，液滴體積逐漸減小直至完全蒸發。液滴在液體基底上鋪展，達到最大鋪展半徑，之后由于蒸發鋪展半徑逐漸減小，蒸發液滴在不混溶液體表面上的鋪展是非釘扎的。在“鋪展主導”狀態下，接觸線的前進速率大于接觸線的后退速率，液滴的直徑隨時間而增大，而液滴體積由于蒸發而減少。這種趨勢導致接觸角減小，接觸線運動阻力增大，從而導致“鋪展主導”狀態向“蒸發主導”狀態過渡。

圖5 h1最小值、h2最大值以及液滴半徑xmax隨時間的演化Fig.5 Temporal evolution of minimal h1，maximal h2，and the drop half-width xmax

2.2 蒸發數的影響

蒸發數E 對鋪展過程的影響如圖6 所示。從圖中可以看出，隨著蒸發數E 的增加，液-氣界面高度逐漸降低，液-液界面高度逐漸升高。蒸發數E表示鋪展時間尺度與蒸發時間尺度的比值，而蒸發時間尺度為液滴自初始靜止時刻開始蒸發至完全干涸所需時間。E 越大，液滴蒸發得越快，質量損失也就越快，相同時刻液滴體積越小，液滴壽命也就越短。驅動液滴鋪展的主要因素包括溫度梯度導致的Marangoni 效應以及蒸發誘導的外向流動。由于蒸發對液滴變形的影響，從液-氣界面蒸發去除的液滴分子必須得到補充以滿足界面力平衡，液滴內部產生從下向上的質量遷移，導致液-液界面不斷升高。

圖6 t = 100時E對液滴形狀的影響Fig.6 The effect of E on the drop shape at t = 100

2.3 表面張力比的影響

圖7為t = 100時不同張力比σ對液滴形狀的影響。如圖所示，隨著張力比的減小，液滴上表面比下表面變形更大[29]，接觸角θ1逐漸減小，θ2逐漸增加。因此，在較小的σ值下，液滴大部分伸出液體基底，而對于較大的σ，液滴則“下沉”到液體基底某一深度。液滴在一定程度上穿透了液體基底，深度取決于表面張力比。

圖7 t = 100時表面張力比對液滴形狀的影響Fig.7 The effect of the surface tension ratio σ on the drop shape at t = 100

2.4 密度比的影響

密度比ρ 的影響如圖8 所示。當液滴密度小于基底液體密度時，液體基底整體向上拱起呈凸形，三相接觸點位于液體基底最高點。而對于ρ ＞1，液體基底整體向下凹陷，三相接觸點位于液體基底的局部最低點，液滴在重力作用下會下沉到液體基底中，密度比越大，下沉深度越大。上層液滴對液體基底的擠壓，導致液體基底厚度從液滴邊緣開始向外逐漸增加。對于負鋪展系數的液滴，考慮到液滴的尺寸足夠小，無論其密度如何，液滴都能夠漂浮[30]。當ρ = 1 情況下，接觸線兩邊基底高度保持不變。顯然，重力是導致液滴呈現圖中所示現象的主要原因。

圖8 t = 100時密度比對液滴形狀的影響Fig.8 The effect of the density ratio ρ on the drop shape at t = 100

2.5 黏度比的影響

黏度比μ 對液滴動力學的影響如圖9 所示。假設液滴上方的氣體具有足夠小的黏度和密度，不會對系統產生影響。圖9(a)表明，液滴中心處的厚度變化過程均要經歷先下降后上升然后再下降三個階段，這是因為重力和毛細力、熱毛細力和蒸發分別在不同階段起主要作用。液滴在蒸發過程中逐漸變薄并產生不規則振蕩運動[6]。早期階段液滴在低黏度比下的流動性更強，導致液滴鋪展更加迅速，對于較大的μ值，黏度比的增加會導致鋪展和收縮速率的降低，同一時刻液滴最大厚度也相對較大。從圖9(b)可以看出，鋪展半徑隨著時間先增加后減小，系統從“鋪展主導”階段逐漸向“蒸發主導”階段過渡。接觸線以恒定速度運動的區域有兩個：一個處于前進階段，另一個處于后退階段[9]。當液滴的黏度增加時，由于黏性耗散和摩擦而產生的阻尼效應會變大[31]，導致鋪展和收縮速率的降低。同時由于蒸發的存在，阻礙了液滴的鋪展，從而縮短了鋪展過程。接觸角的時間演化取決于液滴的鋪展和蒸發以及液滴內部的質量傳輸。對于低黏度液體，初始運動主要由表面張力和初始液滴形狀影響，由于黏度會導致表面張力的變化，因此液滴的初始接觸角是不同的。圖9(c)表明在鋪展的早期階段，隨著液滴鋪展，接觸角減小，當減小至某一極限時，液滴前緣變得不穩定。結果表明，兩種液體間的黏度比對最大液滴厚度、鋪展半徑和接觸角均有一定影響。模擬所得接觸角隨時間的變化趨勢與Sun 等[9]將體積10 μl 甲苯透鏡放置在21℃去離子水表面后的觀察結果相同，均為鋪展階段接觸角減小較快，在蒸發過程中大部分時間保持不變，蒸發過程遵循恒定的接觸角模式。

圖9 黏度比對最大液滴厚度(h2 - h1)max、鋪展半徑xmax和接觸角θ1和θ2的影響Fig.9 The effect of the viscosity ratio on the maximal drop thickness(h2 - h1)max，spreading radius xmax and contact angles θ1 and θ2

2.6 界面溫度分布

通過潤滑理論忽略量級較小的項，消去能量方程中的非穩態項，并計算液滴溫度隨高度的相對變化率。由式(18)和式(19)可知液-液界面(z=h1)和液-氣界面(z=h2)溫度分布和蒸發量為：

式(29)表明，蒸發通過影響液滴界面的溫度分布進而影響表面張力以及液滴鋪展。無論液滴是浸沒還是漂浮在液體基底上，在垂直于基底的溫度梯度下溫度場需關于z 軸對稱。空氣中的熱通量和動量通量較小，對液滴和液體基底的溫度場和速度場無顯著影響，液滴內部溫度波動較小。Buffone[12]進行了乙醇表面FC-72 液滴的蒸發鋪展實驗，并用紅外熱像儀測得了FC-72 液滴表面溫度的分布圖。為了與Buffone[12]的實驗結果進行對比，根據其文獻中提供的物性參數計算了密度比ρ=2.1，熱導率比λ=0.34，黏度比μ=0.43，表面張力比σ=0.46 時，液滴蒸發演化過程中的液-氣界面溫度分布和液-液界面溫度分布以及溫度梯度分布，如圖10所示。圖10(a)為不同時刻液滴界面溫度分布，由于溫度與高度呈反比，故液-液界面溫度高于液-氣界面，且在液滴中心處達到峰值，液-氣界面則與之相反。圖10(b)、(c)表明液-液界面的溫度梯度大于液-氣界面。由于邊緣處蒸發更快，蒸發冷卻效應使得液滴邊緣較中心區域溫度更低。計算所得液-氣界面溫度的數量級與Buffone[12]實驗結果無量綱化之后液滴表面溫度數量級一致，且液滴表面溫度隨時間的變化趨勢也與實驗結果一致，即表面溫度均隨時間而升高。

圖10 液-液和液-氣界面溫度與溫度梯度的演化過程Fig.10 Evolutions of liquid-liquid and liquid-gas interfaces temperatures and temperature gradients

2.7 與實驗結果的比較

Dussaud 等[11]在23℃時觀測了表面張力比分別為1和0.4，密度比分別為1和0.9的辛烷和甲苯液滴在去離子水表面的蒸發行為，發現液滴前緣以t1/2的冪律向前推進，即xmax- t1/2。為了與實驗進行對比，采用相同的表面張力比和密度比計算所得的不混溶液體表面上蒸發液滴鋪展早期階段鋪展半徑xmax隨時間變化關系如圖11 所示(箭頭指向σ 增大的方向，圖中的實驗結果采用本文的無量綱定義式進行了無量綱化)。圖11 表明，實驗與數值模擬結果具有相同的變化趨勢，但存在一定的誤差，原因可能在于模擬液滴是二維的，無法完全與實驗進行比較。此外，由于實驗未涉及加熱，而加熱會使液滴蒸發加快進而導致鋪展速率減小，故在鋪展的早期階段，鋪展速率模擬值比實驗值要小，隨著時間的進行，蒸發對液滴鋪展的影響越來越大，在鋪展的后期階段鋪展半徑模擬值小于實驗值。

圖11 不同表面張力比σ時鋪展半徑的時間演化對數圖Fig.11 Log-log plot of the temporal evolution of spreading radius for different surface tension ratio σ

3 結 論

采用潤滑理論推導了液-氣和液-液界面的演化方程。 通過數值模擬研究了上層液滴外形以及鋪展半徑的時間演化規律，以及蒸發、重力、毛細力和Marangoni效應之間的耦合作用對接觸角、鋪展半徑、溫度場和界面幾何形狀的影響。模擬結果表明演化過程存在兩種狀態，一種是“鋪展主導”的狀態，另一種是“蒸發主導”的狀態。隨著時間的推移，接觸角逐漸減小，接觸線運動阻力增大，“鋪展主導”狀態逐漸向“蒸發主導”狀態過渡。蒸發量與液滴高度呈反比，蒸發數越大，液滴蒸發得越快。液滴高度隨時間逐漸降低，表面溫度則越來越高。液滴在一定程度上會伸出和穿透液體基底，深度取決于表面張力比和密度比。低黏度比下的流動性增強會導致更快的鋪展，黏度比的增加會導致鋪展和收縮速率的降低。液滴前緣首先沿液體基底表面前進，在達到最大接觸直徑后由于蒸發而逐漸后退。在鋪展的早期階段，鋪展半徑xmax隨時間變化關系與實驗中得到的xmax- t1/2規律具有相同趨勢。