豎直平板上含活性劑液膜排液過程不穩(wěn)定現(xiàn)象的數(shù)值模擬

江 蕊，李春曦，劉成治，葉學民

（河北省低碳高效發(fā)電技術重點實驗室（華北電力大學），河北 保定 071003）

0 引言

浸涂是將基板從溶液池中浸漬提拉，經(jīng)排液干燥，使表面形成一層均勻薄膜的過程[1]。該工藝的作用是保護基板材料、改善表面性能，故而被廣泛用于工業(yè)領域[2,3]。

在排液過程中，液膜表面的均勻性和穩(wěn)定性是影響涂層質量的關鍵因素[4]。深入研究該過程中液膜的流動特性和不穩(wěn)定現(xiàn)象的發(fā)展機理，對提高涂層質量具有重要意義。

目前，關于平板液膜排液不穩(wěn)定現(xiàn)象的理論研究大多采用線性穩(wěn)定性分析[5,6]，確定出液膜的臨界穩(wěn)定狀態(tài)及最危險波數(shù)，而很少直接模擬其不穩(wěn)定現(xiàn)象。雖然穩(wěn)定性分析能篩選出令液膜失穩(wěn)的因素，但為驗證理論分析結果，還需通過數(shù)值模擬再現(xiàn)平板液膜排液中的不穩(wěn)定現(xiàn)象。

與平板液膜排液具有相似數(shù)理模型的線框液膜排液問題已有成功算例可供借鑒。文獻[7]對附著于線框表面的液膜建立了二維排液模型，在初始條件中加入余弦波展向擾動，模擬了液膜表面的不穩(wěn)定現(xiàn)象，發(fā)現(xiàn)Marangoni效應和表面黏性會影響液膜流動的穩(wěn)定性。文獻[8]針對線框液膜排液的不穩(wěn)定現(xiàn)象進行了進一步分析，指出：活性劑分布的不均勻性隨排液發(fā)展而加劇，Marangoni效應對擾動發(fā)展具有重要影響。

上述研究表明：Marangoni效應與液膜表面的不穩(wěn)定現(xiàn)象緊密相關。雖然平板液膜的排液模型與線框具有相似性；但相比線框排液過程，平板液膜還需滿足表面無滑移和無通量等邊界條件[9]。此外，對于含活性劑液膜，其排液過程中可能出現(xiàn)的不穩(wěn)定性現(xiàn)象及其內在機理也尚待研究。

為此，本文擬建立平板表面含活性劑液膜的排液模型，推導表面速度、液膜厚度及活性劑濃度的二維方程組；通過在彎月面處引入隨機波初始擾動來模擬液膜表面的不穩(wěn)定現(xiàn)象，分析液膜流動特性、Marangoni效應對擾動過程的影響以及不穩(wěn)定現(xiàn)象的發(fā)展機理，以期為提高浸涂工藝中液膜的涂覆質量提供理論基礎。

1 理論模型建立

自液池中提拉出的豎直平板表面會附著一層液膜，液膜受重力作用排液，過程如圖1所示。圖1中，平板懸垂靜止，液膜頂端厚度為零，底端與液池相連形成彎月面；x、y、z分別為展向、流向及與平板垂直的方向。假設液膜表面含有不溶性活性劑，假設流體為不可壓縮黏性流體，并忽略排液中液膜蒸發(fā)的影響。排液過程中，平板左右兩側的液膜對稱，因而僅需分析其一側液膜的變化特征。

圖1 豎直平板上液膜排液過程Fig. 1 Film draining on a vertical plate

1.1 控制方程

排液過程的控制方程包括連續(xù)方程和運動方程：

式中：含“*”的物理量表示有量綱值；“·”表示2個向量的內積； ?*為哈密頓算子； ?*2為拉普拉斯算子；V*=(u*,v*,w*)；j=(0,1,0)；μ*為動力黏度；ρ*為密度；p*為壓力；t*為時間；g表示重力加速度。

在z*=0處，邊界條件為：

在z*=h*處，運動學條件為：

液膜氣液界面滿足應力平衡邊界條件[7]：

式中：“×”表示2個向量相乘的外積；“：”表示張量的并聯(lián)雙點積；n*、t*分別代表液膜表面的單位外法向和切向矢量；T*為液膜應力張量；H*為曲面平均曲率；σ*為表面張力；κs*、μs*分別為液膜表面的膨脹黏度和剪切黏度為液膜的法向曲率；為表面單位張量。

假設表面張力與活性劑濃度滿足線性關系[10]：

在液膜自由面z*=h*(x*,y*,t*)處，活性劑濃度的輸運方程為[11]：

1.2 方程無量綱化

方程無量綱變換采用式（9），變換后得到的無量綱物理量不帶上標“*”。

式中：ε=d*/l*；d*為液膜特征厚度；l*為動態(tài)彎月面長度，l*=(d*D*2)1/3，其中D*為彎月面半徑；特征速度U*=ρ*gd*2/μ。S為液膜表面黏度；S1為剪切黏度；M為Marangoni數(shù)；Pe為Peclet數(shù)，表征活性劑擴散程度。

將式（9）代入式（1）～（8）中進行無量綱化。結合潤滑理論（ε<<1），忽略數(shù)量級為ε2的小量，可得無量綱后x和y方向上表面速度、液膜厚度和活性劑濃度演化方程為：

2 數(shù)值模擬

本模型的求解區(qū)域設定為Lx=12π，Ly=12π；計算網(wǎng)格劃分為100×100；時間步長設定為0.01。所取參數(shù)和范圍參照了本課題組前期文獻[11]，如表1所示。

表1 典型尺度參數(shù)和無量綱參數(shù)取值或變化范圍Tab. 1 Value or variation range of typical scale parameters and dimensionless parameters

采用“FreeFEM++v4.7”求解方程組（10）～（13）。該軟件在高精度數(shù)值求解液膜流動問題中具有顯著優(yōu)勢[12,13]。

計算中邊界條件設置如下：

平板頂端y=0處，厚度為0，無液體流動，即：

底端y=Ly處，液膜與液池相接，形成彎月面。根據(jù) Young-Laplace方程[14]，液膜厚度的一階和二階導數(shù)滿足以下關系：

靠近液池區(qū)域的液膜沿y方向無活性劑流動，即：

平板展向兩側采用相同的邊界條件：

含活性劑液膜排液過程的數(shù)學描述是一個與初始條件相關的初值問題，將擾動施加在液膜厚度和濃度的初始條件中[15]。考慮到實驗過程中擾動往往隨機產生，并未出現(xiàn)確定的主導波數(shù)，故本文參考了文獻[16]中描述含活性劑液膜水平鋪展過程中不穩(wěn)定現(xiàn)象所采用的擾動形式。本研究中，為模擬豎直平板上排液過程中液膜表面出現(xiàn)的展向擾動波，在液膜厚度初始條件的彎月面附近yd=35處添加隨機波數(shù)擾動，如式（18）所示：

活性劑濃度的初始條件如式（19）：

式中：Rand[x,y]表示[x,y]閉區(qū)間內的隨機數(shù)；yd是彎月面位置，取yd=35；y1=36.69；y1

式（19）的含義為：液膜主體部分活性劑濃度為0，自接近彎月面處（y≈33）活性劑濃度急劇增加到1，與液池中濃度相同，且在yd=35處引入隨機擾動。

在式（19）中，人為施加了一極大的活性劑濃度梯度，其結果與實際工程中快速提拉平板形成的液膜表面活性劑的濃度分布接近。

3 結果與分析

彎月面附近為液膜最薄區(qū)域，當其受到微小擾動后極易出現(xiàn)不穩(wěn)定特征。為探究局部擾動能否隨排液持續(xù)而維持不穩(wěn)定性的問題，在初始條件中施加展向擾動后：需要首先分析液膜厚度和流動特性的變化，以此驗證平板液膜二維模型的合理性；然后模擬液膜表面的不穩(wěn)定特征，揭示局部展向擾動的成因及發(fā)展過程。

3.1 M=1時液膜流動特征

當M=1、S=10、S1=1時，模擬所得x=15處的液膜厚度和流向表面速度如圖2所示，液膜厚度和活性劑濃度的二維演化過程如圖3所示。色條顏色反映了活性劑濃度高低。

圖2 M=1時主流方向的液膜厚度及表面速度分布Fig. 2 Film thickness and surface velocity distribution in the mainstream direction when M=1

結合圖2（a）和圖3知，重力作用促使液膜向下排液，上端液膜厚度逐漸減少，彎月面以下厚度則快速增加。隨排液進行，動態(tài)彎月面位置上移到y(tǒng)=32.5處，液膜厚度變薄、速率減緩；自t=80開始，厚度幾乎不變。在圖2（b）中，當t=20時排液速度較初始時刻大幅降低，底端彎月面附近的表面速度首先出現(xiàn)負值，這是因初始條件中的大濃度梯度誘發(fā)了 Marangoni效應，從而使液膜逆流；t=40時逆流范圍有所擴大，逆流速度提高，同時逆流也攜帶底部的活性劑向上輸運，減弱了流向的Marangoni效應，因而后續(xù)幾個時刻的表面速度呈現(xiàn)出主流速度和逆流速度同時減小的特點。圖 2所示排液過程中，液膜厚度變化及逆流特征與本課題組研究一維平板浸涂過程所得的結果[9]4-6一致。

圖3 二維液膜表面的厚度和濃度分布Fig. 3 Thickness and concentration distribution of 2D film surface

由于本文重點關注彎月面附近區(qū)域（y值位于30～35）液膜表面的不穩(wěn)定現(xiàn)象，故截取液膜厚度0

由圖3可知，在排液過程中，液膜沿流向變化的同時也沿展向變化。在t=0初始時刻，活性劑集中在液膜底部，在彎月面位置處存在著具有隨機擾動波幅大小的擾動波紋；t=20時，逆流促使活性劑分布較均勻，擾動波紋消退，彎月面附近沿展向相對均勻；自t=60時開始，活性劑濃度出現(xiàn)局部峰值，展向的濃度分布不再均勻，彎月面處衍生出微弱的起伏波紋；而到排液末期t=100時，活性劑濃度分布的不均勻度更加顯著，彎月面處波紋起伏程度也明顯加劇，形成多個凸起結構。

3.2 液膜表面的不穩(wěn)定特征

為進一步考察液膜表面不穩(wěn)定現(xiàn)象的演化過程，參考光學干涉法測量液膜厚度的實驗原理，將液膜厚度在平板上自上而下的變化用漸變色條表示，色條沿流向的疏密程度反映液膜厚度大小，沿展向的彎曲度表示擾動發(fā)展特征，結果如圖 4所示。

圖4 M=1時液膜表面擾動發(fā)展過程Fig. 4 Disturbance evolution of the film surface when M=1

由圖4知，在t=0時，與初始條件設置的隨機波數(shù)擾動相對應，液膜底端彎月面處出現(xiàn)了擾動波紋；但在t=20時，擾動波紋消退，液膜表面在展向上呈現(xiàn)相對平整的狀態(tài)。排液后期自t=60時開始，液膜表面在y≈27處出現(xiàn)微弱波紋，且隨排液進行擾動波紋不斷發(fā)展；到t=100時，條紋紊亂，擾動波紋起伏劇烈，在y=30、x位于20～30的區(qū)域出現(xiàn)了2個液膜厚度的高峰。圖中液膜表面的條紋特征反映了排液過程中擾動波紋的發(fā)展動態(tài)。

在排液初期與后期，液膜表面的擾動特征具有明顯的差異性，這表現(xiàn)在：t=0排液初期，在厚度初始條件中施加擾動波紋的波幅大小隨機，其影響范圍僅限于底部彎月面y>35位置，且擾動迅速消退；而在t>60排液后期，重新出現(xiàn)的擾動波紋波幅較小，但卻覆蓋了y值在15～30的大范圍區(qū)域；隨著對排液進行擾動的波紋波幅的增大，擾動的影響范圍不斷擴展；到t=100時，存在擾動的區(qū)域基本覆蓋了y值位于10～35的絕大部分液膜表面。

局部擾動在排液后期液膜極薄的狀態(tài)下發(fā)展更為充分，并隨排液逐漸演化為全局擾動。

3.3 不穩(wěn)定現(xiàn)象分析

3.3.1 Marangoni效應的影響

在S=10、S1=1時，M分別取為1、10和30，分析 Marangoni效應對液膜排液不穩(wěn)定現(xiàn)象的影響。

圖5示出了不同M值時液膜表面速度演化過程。由圖 5（a）所示的x=15處的流向速度vs分布可知：vs分布在初始時刻并無差別；但在t=60時，出現(xiàn)大小不同的逆流區(qū)；逆流區(qū)隨M增大而擴大（如M=30時，對應的逆流區(qū)上邊界接近y=20；液膜逆流區(qū)占據(jù)近一半的流動域），逆流速度的最大值也隨M增大而提高；在t=100排液末期時，Marangoni效應相對t=60時減弱，導致除彎月面區(qū)域外，M=1和M=10所對應的逆流區(qū)消失，M=30所對應的逆流區(qū)也縮小至 15

圖5 不同M時液膜表面速度演化過程Fig. 5 Evolution of film surface velocity at different M

圖5（b）所示為y=35處液膜表面展向上的速度分布。

（1）在t=0時，展向速度us為0，液膜受重力主導沿著豎直方向流動。

（2）當t=60時，液膜表面衍生出展向速度。M=1情形下，us≈0.03，且沿展向相對均勻。隨M增大，us均勻性變差。當M=30時，us在x值位于10～25范圍內出現(xiàn)多個峰谷，且速度在正負數(shù)值間變化——這表明沿展向上的活性劑濃度分布不同，液膜受到了不同方向的 Marangoni作用。此時，液膜沿流向的逆流現(xiàn)象依舊顯著，見圖5（a）；而展向上，因 Marangoni效應，液膜表面也呈現(xiàn)出速度不均勻的特征。

（3）在t=100時，us出現(xiàn)震蕩特征，不同x位置處的速度數(shù)值差異顯著。隨M增大，us峰值和均值呈增大趨勢，液膜呈現(xiàn)如圖4所示的波浪起伏特征。

排液后期，液膜表面逐漸發(fā)展的擾動波紋并非只出現(xiàn)在底端彎月面附近，而是覆蓋了絕大部分流動域（y>20）。當t=100時，液膜表面的擾動現(xiàn)象尤其明顯，如圖6所示。

圖6 t=100時刻液膜表面不穩(wěn)定現(xiàn)象Fig. 6 Instability of film surface at t=100

另外，當M值較大時，液膜厚度在流向和展向上均發(fā)生較大的起伏變化，擾動波紋尺度增大，不規(guī)則性增強。特別是在M=30時，波紋起伏最為突出，這表明液膜處于極不穩(wěn)定的流動狀態(tài)。

3.3.2 擾動現(xiàn)象的成因

排液過程初始時刻，彎月面附近區(qū)域因施加擾動而出現(xiàn)不穩(wěn)定波紋，液膜受重力作用豎直向下流動，擾動的展向發(fā)展受到阻礙，因而擾動發(fā)展微弱甚至消退。隨排液進行，活性劑逆流使液膜底端（y>33）處的活性劑濃度相對初期減小，而中上部濃度從零逐漸增加（見圖7），因而彎月面附近區(qū)域（30

圖7 M=30液膜活性劑濃度分布Fig. 7 Distribution of film surfactant concentration when M=30

為進一步量化不穩(wěn)定現(xiàn)象，定義瞬時擾動能量He如式（20）：

式中：hb為未施加擾動時的基態(tài)液膜厚度。

某一時刻液膜表面擾動總能量Eh如式（21）：

圖8示出了液膜表面擾動能量Eh隨時間的演化過程。

圖8 液膜表面擾動能量隨時間變化Fig. 8 Variation of film disturbing energy with time

（1）t<20排液前期，人為施加的擾動在彎月面附近（30

（2）到t>60排液后期，彎月面h<0.4h0處，受 Marangoni力作用，彎月面區(qū)域整體活性劑濃度低，濃度中施加的擾動影響強烈。該區(qū)域擾動波紋隨排液進行擴展到大部分液膜區(qū)域。整個液膜表面擾動能量大幅增加，如Et=100≈20Et=0。

排液后期，彎月面區(qū)域液膜厚度小，且存在低濃度的活性劑分布，施加的展向擾動影響突出。由圖8知，隨M增大，t>60后液膜擾動能量增長的速率加快；尤其是在M=30時，擾動能量增長速率較M=1時大幅增加。這是由于當M增大時，活性劑沿流向的逆流區(qū)域擴大；又因展向擾動作用，活性劑濃度梯度誘發(fā)的 Marangoni效應也強烈作用于展向方向，于是：擾動波紋的影響不局限于彎月面位置，而是向整個液膜表面發(fā)展，從而導致液膜表面不穩(wěn)定性增強，整體擾動能量大幅增長。

4 結論

（1）平板液膜排液的特點為：排液前期受重力主導液膜厚度迅速減薄，平板底端活性劑濃度梯度誘發(fā)的 Marangoni效應促進液膜逆流，流向速度出現(xiàn)負值；后期活性劑濃度分布呈現(xiàn)二維不均勻性，液膜表面衍生出展向速度，沿展向出現(xiàn)多個峰谷，液膜彎月面區(qū)域出現(xiàn)褶皺的凸起結構。

（2）施加局部擾動后，液膜表面擾動發(fā)展的機制為：初始時期，彎月面區(qū)域因施加了擾動條件而出現(xiàn)不穩(wěn)定波紋，但因重力流動的阻礙作用，擾動的展向發(fā)展受限以至于快速減弱；排液后期，彎月面液膜厚度極薄，該區(qū)域內活性劑濃度整體較低，流向和展向方向均存在 Marangoni效應，導致擾動波紋再現(xiàn)并向頂端擴展。

（3）施加展向擾動的排液過程中，尤其是在排液后期液膜厚度極薄的狀態(tài)下，Marangoni效應強化了液膜表面的擾動特征。隨 Marangoni效應增強，液膜表面擾動波紋擴展的范圍擴大，波紋起伏不規(guī)則性增強，液膜表面整體擾動能量大幅增長，這對液膜的穩(wěn)定性產生不利影響。