邵明玉,張忠偉,馬馳騁,荊 棟,華 珍

(山東理工大學 交通與車輛工程學院,山東 淄博 255000)

薄膜流動現(xiàn)象普遍存在于化學工程、生物工程、機械潤滑等領域,作為一類典型的流體動力學問題,蘊含著豐富的非線性行為,引起了廣大研究者的關注[1-3]。當薄膜在平面、圓柱面或者沿著細絲流動時,含運動接觸線的薄膜流動存在多種有趣的形態(tài),比如穩(wěn)定行波演化,或指突不穩(wěn)定。影響薄膜流動特性的因素很多,如壁面拓撲結構[4]、靜電場力[5]、熱毛細力[6]和滑移邊界效應[7]等,而且在薄膜流動過程中往往多種因素共同作用。

從流體動力學角度來講,薄膜流動屬于低速開式的交界面流動,自由面處的非線性動力學行為分析和數(shù)值模擬一直是薄膜流動研究的重點之一。相對于理論分析和實驗驗證,數(shù)值模擬在薄膜下落問題研究中更加普遍,而且數(shù)值方法眾多。胡軍等[8]采用拉格朗日有限元法,開展了均勻加熱下落薄膜的直接數(shù)值模擬；Lin等利用有限差分方法離散求解空間,用alternating direction implicit(ADI)方法進行時間迭代,雖然精度很高,但算法施加相對繁瑣[9-10]；潘曉軍等利用傅里葉偽譜法與Runge-Kutta法相結合開展了水平基底上薄膜流體流動的數(shù)值模擬[11-12]；Gao等基于有限體積法研究了薄膜流動[13]；Gu等采用計算流體力學方法研究了沿斜板下落薄膜的流動規(guī)律,并分析了兩相流工況下的薄膜流動形態(tài)[14]。

近年來,更加新穎的仿真技術手段應用到流體動力學仿真中。李杰和Raach等采用OpenFoam軟件實現(xiàn)了薄膜流動的二維數(shù)值模擬,獲得了大量的數(shù)值結果[15-16]；Wehinger利用STAR-CCM+模擬軟件開展了液膜下降流動的數(shù)值模擬,通過局部網(wǎng)格加密獲得了比較精確的仿真結果[17]；Hu使用有限單元求解器Dolfin開展了非牛頓流體的三維薄膜流動仿真,并分析了流動的穩(wěn)定性[18]。從文獻檢索中可以發(fā)現(xiàn),盡管有限單元法早在1984年便開始用于下降薄膜的數(shù)值模擬,但是在流體動力學仿真中,有限差分方法和有限體積法的使用更加廣泛。隨著薄膜流動問題的深入研究,尋求一種更為簡單有效的數(shù)值方法將會促進該類問題的研究。

本文在長波近似的理論框架下,推導了倒置平面基底上的薄膜流動潤滑方程,然后針對該方程建立有限元模型,使用開源有限元軟件Freefem++仿真薄膜二維流動和三維流動中的運動接觸線演化過程,基于流動形貌的演化過程,系統(tǒng)地分析了倒置平面基底上的液膜流動演化規(guī)律。

1 薄膜運動方程

1.1 基本方程

3種沿平板下落的薄膜流動現(xiàn)象如圖1所示,分別為沿斜板流動、豎直流動和倒置基底流動。這3種流動都屬于重力驅(qū)動流動,前兩種現(xiàn)象已經(jīng)得到了非常廣泛的研究,在前人工作的基礎上,本文重點研究倒置平面基底上的薄膜流動問題。研究中假設流體為黏性不可壓縮的牛頓流體,因此，滿足連續(xù)方程和Navier-Stokes方程。

圖1 平面薄膜流動Fig.1 Falling films on planes

(1)

(2)

(3)

(4)

圖1中所示的薄膜與基底的接觸面采用無滑移邊界,即

u=v=w=0。

(5)

在自由界面處滿足動力方程和張力平衡方程，

(6)

(7)

1.2 潤滑理論

基于薄膜流動特征,采用潤滑理論推導其流動方程。由于薄膜厚度相比其他兩個方向上尺寸較小,因此，引入小參數(shù)ε=h0/L。其中：h0為薄膜平均厚度；L為其他方向的特征長度。無量綱坐標可以表示為

(8)

令U表示x方向的速度尺度,速度項的無量綱表達式為

(9)

通過特征長度和特征速度,可以確定特征時間,同時將體積力勢函數(shù)和壓強定義為

(10)

(11)

將無量綱參數(shù)表達式(8)～(11)帶入連續(xù)方程和Navier-Stokes可以得到

(12)

(13)

(14)

(15)

其中，Re表示雷諾數(shù)。為簡化表達，下文省略無量綱符號。本文主要研究重力驅(qū)動的薄膜流動,體積力只考慮重力(ρgsinα,0,-ρgcosα),略去式(13)～式(15)中包含ε及其高階小量的部分，可以得到

(16)

(17)

(18)

對式(18)在z方向上積分,代入邊界條件式(7)可以得到壓強函數(shù)為

p=ρgcosα(h-z)-σκ+p0。

(19)

其中，平均曲率為

κ≈2h。

(20)

得到壓強的表達式(19)后,對式(16)和式(17)在z方向上積分兩次,利用邊界條件可以得到x方向和y方向的速度表達式，

(21)

(22)

動力邊界條件(6)也可以用積分形式表示，

(23)

通過平均速度進一步簡化式(21)和式(22)可得

(24)

(25)

將式(24)和式(25)代入式(23)，可得倒置平面基底上薄膜流動的潤滑方程，

(26)

根據(jù)表面張力和重力項的平衡條件[19],液膜厚度的演化控制方程可以進一步化簡為

(27)

其中，Dθ=(3Ca)1/3cotθ,Ca=μU/σ代表毛細管數(shù),式(27)雖然比式(26)簡單,但仍然是強非線性、高階偏微分方程,采用差分方法求解該方程依然會比較復雜。

2 有限元方程

根據(jù)有限單元法和Freefem++編譯環(huán)境[20],求解微分方程時，需要首先建立微分方程的變分形式。直接建立方程(27)的變分形式比較困難,為了得到方程(27)的微分形式,引入一個新的參數(shù)φ，

φ=2h。

(28)

式(27)可以改寫為

(29)

式(27)的求解轉(zhuǎn)換為式(28)和式(29)的求解,雖然方程數(shù)目增加,但是求解效率更高。通過分部積分,推導出式(28)和式(29)的變分形式為

(30)

(31)

其中：H和ψ分別為h和φ的試函數(shù)；S表示計算域。

初始條件和計算域如圖2所示,在Freefem++中對y=0和y=Ly處施加周期邊界條件,初始條件為光滑函數(shù)連接的兩部分薄膜厚度，

(32)

其中：b為前驅(qū)膜的厚度；x0為y的函數(shù)，

x0=xf-A0cos(2πy/λ)。

(33)

其中：A0和λ分別為初始擾動的幅值和波長；xf表示連接處位置。

3 倒置基底薄膜流動二維模擬

在Freefem++中無法直接進行二維仿真,但是可以通過設置一個較小的y方向尺寸近似開展二維仿真,忽略y方向上函數(shù)值的變化，得到的結果與二維模擬十分接近,完全能滿足計算精度。有限元模型和邊界條件如下圖2所示。使用前驅(qū)膜假設,假定固體表面完全潤濕,前驅(qū)膜厚度為b,液膜從左側計算域流進,在x=0處厚度設置為1,液膜初始形貌根據(jù)式(32)給定,仿真模擬中xf為15。

圖2 有限元網(wǎng)格和邊界條件Fig.2 Finite element grids and boundary conditions

倒置平面基底上的液膜流動形貌由重力和表面張力的共同作用決定。傾斜角度越大,作為驅(qū)動力的重力作用越強,相較于表面張力的作用更大,從而更易引起液氣交界面的非線性演化。文中Dθ取值為-0.61,-1.56和-1.93, 分別代表傾斜角度122°, 148°和153°。 首先， 計算Dθ=-0.61時的液膜流動鋪展過程,有限元模型和不同時刻液膜的鋪展演化形貌如圖3所示,可以看出,液膜厚度在y方向上的變化完全可以忽略,說明采用這種近似方式開展液膜流動的二維仿真是可行的。

圖3 二維問題的近似求解域和數(shù)值演化結果Fig.3 Approximate calculation domain for two dimensional problems and time evolution results

只保留圖3中y=0.5處的數(shù)值結果,即得到如圖4中所示的二維演化形狀。從圖4中可以看出,當Dθ=-0.61,傾斜角度較小時,液膜流動呈現(xiàn)出一種穩(wěn)定的行波態(tài),經(jīng)過一段時間的演化,毛細脊的高度不再發(fā)生改變,穩(wěn)定地向前推進。可以推斷出,當?shù)怪媒嵌容^小時,液氣接觸表面的鋪展是運動穩(wěn)定的,其流動穩(wěn)定性與豎直基底上的液膜流動穩(wěn)定性一致。當Dθ=-1.56,液膜流動規(guī)律變得十分復雜,t=10～30時間內(nèi),前進波呈現(xiàn)出一種強烈的衰減形式,隨著時間增長,液膜呈現(xiàn)出一種強衰減形態(tài)和孤立波并存的形式,靠近接觸線前緣逐漸演化出一系列孤立波,而靠近左側邊界處依然存在強衰減的前進波,如圖4B中t=40～60時所示。在液膜鋪展過程中,液膜厚的區(qū)域遷移快,而液膜薄的區(qū)域遷移慢,因此，在移動接觸前緣區(qū)域能夠觀察到兩個相鄰的波峰融合成一個大的孤立波形態(tài)。而當Dθ=-1.93時,液膜運動的不穩(wěn)定性更強,在較短的時間內(nèi)就發(fā)展成多個孤立波存在的形式,如圖4C所示。通過二維的薄膜流動仿真模擬可以確定,θ越大,液膜的流動越不穩(wěn)定,主要原因是θ越大，重力驅(qū)動作用越強,導致表面張力無法約束液氣接觸面的穩(wěn)定鋪展。

圖4 不同倒置角度基底薄膜流動二維演化Fig.4 Two dimensional flow simulation of the falling film on an inverted plane for different angles

4 倒置基底薄膜流動三維模擬

在液膜流動過程中,液氣兩相界面往往會發(fā)生失穩(wěn)以指狀形貌流動鋪展,二維仿真可以觀察到兩相界面的起伏，但是無法研究其黏性指進行為。因此,進一步針對方程(27)開展三維有限元數(shù)值模擬。三維演化模擬結果如圖5～圖7所示。計算域為(0≤x≤100,0≤y≤40)。為了考察液膜流動的三維演化特性,施加初始條件為

x0=xf-Aicos(2πy/λi)。

(34)

其中：i=1，2，…，50；Ai和λi分別表示第i階的擾動幅值和擾動波長。

圖5～圖7分別給出了Dθ為-0.61,-1.56和-1.93這3種工況下, 倒置平面基底上薄膜流動的三維演化過程。 對比圖4和圖5～7可以看出, 液膜的三維流動和二維模擬具有密切的聯(lián)系性, 二維流動反映了三維流動的一些基本特性, 如移動接觸區(qū)域的遷移速度、 流動形貌等,通過測量發(fā)現(xiàn)， 圖4A和圖5中毛細脊的移動速度和高度是完全一致的, 圖4B和圖4C中二維演化形貌也分別與圖6和圖7的演化形貌在x方向上是相對應的。

圖5 Dθ=-0.61工況下不同時刻下落液膜的三維演化模擬Fig.5 Three dimensional simulation of time evolution of the falling film for case of Dθ=-0.61

圖6 Dθ=-1.56工況下不同時刻下落液膜的三維演化模擬Fig.6 Three dimensional simulation of time evolution of the falling film for case Dθ=-1.56

圖7 Dθ=-1.93工況下不同時刻下落液膜的三維演化模擬Fig.7 Three dimensional simulation of time evolution of the falling film for case Dθ=-1.93

從圖5可以看出,當Dθ=-0.61時,液膜以固定的前進速度、穩(wěn)定的行波形式向下游移動,并沒有發(fā)生黏性指進現(xiàn)象，初始擾動Aicos(2πy/λi)

完全消失, 說明這種流動形式是穩(wěn)定的。如圖6所示，當Dθ=-1.56時,除了能得到x方向的演化形態(tài),也可以得到y(tǒng)方向的演化形態(tài),在t=20和t=40時,自由液面以強衰減形態(tài)為主,而經(jīng)過長時間的演化,在t=60時出現(xiàn)了大量的孤立波,與圖4B中的結果是一一對應的。經(jīng)過一段時間演化,液氣兩相接觸面在y方向上厚度出現(xiàn)了明顯變化,即發(fā)生了黏性指進現(xiàn)象,說明當傾斜角度較大時,自由液面的運動是不穩(wěn)定的。當傾斜角進一步增大,即Dθ=-1.93時,液膜流動的不穩(wěn)定性更強,如圖7所示，可以看出,經(jīng)過較短時間的演化,在t=30時發(fā)生了明顯的黏性指進現(xiàn)象,強衰減的流動模態(tài)基本消失。從圖6和圖7可以發(fā)現(xiàn),孤立波是指形結構的重要幾何特征,隨著時間的增加,孤立波會成為脫落的液滴,導致計算程序的不收斂,這是由于有限單元法對連續(xù)性的要求所決定的。

5 結語

本文通過理論分析并結合數(shù)值模擬,主要研究了倒置平面基底上含運動接觸線的薄膜流動問題。

1) 基于潤滑近似理論,建立了倒置平面基底上薄膜流動的潤滑方程,采用有限元方法基本原理,引入新參數(shù)將薄膜流動演化方程轉(zhuǎn)變?yōu)閮蓚€方程,推導了薄膜流動強非線性高階偏微分方程的變分形式。通過開源有限元軟件Freefem++對薄膜流動運動接觸線的演化開展了二維數(shù)值模擬和三維數(shù)值模擬,它們具有良好的對應性。

2) 數(shù)值結果表明，對于倒置基底上的薄膜流動,倒置角度越大,重力的作用越強,自由液面流動的不穩(wěn)定性增強。傾斜角度較小時,液膜以穩(wěn)定的行波態(tài)運動,當傾斜角度變大,可以觀察到強衰減形態(tài)和孤立波。對比三維鋪展形貌和二維形貌可以發(fā)現(xiàn)，孤立波是黏性指進的重要幾何特征,因此，可以說出現(xiàn)孤立波意味著液氣接觸面有黏性指進現(xiàn)象的發(fā)生。

本文采用的有限元分析步驟對于分析該類問題具有一定的推廣性,由于薄膜流動的運動微分方程都具有類似的形式,因此，容易將本文的工作擴展到其他類似問題中。