袁 浩，黃宏寶，郭平措，何小瀧

(1.重慶交通大學(xué)重慶西南水運(yùn)工程科學(xué)研究所，重慶 400074；2.四川大學(xué)水力學(xué)及山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065；3.青海省水利水電勘測(cè)設(shè)計(jì)研究院，青海 西寧 810001)

液滴撞擊液膜現(xiàn)象廣泛存在于泄洪霧化、噴墨打印、雨滴沖擊地面積水等自然和工程現(xiàn)象中。通常液滴撞擊固壁邊界后在固壁上發(fā)生聚集現(xiàn)象形成液膜，隨后發(fā)生液滴撞擊液膜現(xiàn)象。早期研究表明液滴撞擊液膜后演化過程主要與液滴碰撞速度、流體性質(zhì)、液膜厚度及壁面物理化學(xué)性質(zhì)相關(guān)。

迄今為止，前人對(duì)液滴撞擊液膜過程進(jìn)行了大量試驗(yàn)及數(shù)值模擬研究。試驗(yàn)研究主要通過高速攝像機(jī)捕捉液冠演化過程，并分析各無量綱參數(shù)對(duì)液冠演化過程的影響，主要包括雷諾數(shù)Re、韋伯?dāng)?shù)We和液膜厚度[1-3]。Rioboo等[3]通過試驗(yàn)研究將液滴撞擊液膜后形態(tài)演化分成了融合、反彈、部分反彈、液冠及飛濺5種情況。試驗(yàn)結(jié)果表明液冠半徑與無量綱時(shí)間的平方根成正比[1-2]，液滴撞擊液膜后在液膜中形成的速度不連續(xù)性是液冠形成的原因[1]。然而液膜并非靜止的，通常會(huì)在平面上發(fā)生流動(dòng)。Alghoul等[4]首先研究了液滴撞擊移動(dòng)液膜的過程，隨后不同參數(shù)對(duì)液滴撞擊液膜的影響研究相繼展開，如Castrejón-Pita等[5]將液滴撞擊移動(dòng)液膜的形態(tài)歸納為融合、劇烈飛濺、部分或完全反彈以及液冠這4種形態(tài)，而形態(tài)演化的主要影響參數(shù)是初始液膜運(yùn)動(dòng)速度和碰撞速度；Gao等[6]研究了液冠延伸率和液冠半徑隨時(shí)間的變化，通過試驗(yàn)研究得到了液滴撞擊移動(dòng)液膜時(shí)液冠半徑與時(shí)間的關(guān)系和液冠縱向延伸率變化規(guī)律。

宏觀數(shù)值模擬可以獲得液冠演化及破碎過程中流場(chǎng)及渦量場(chǎng)的詳細(xì)變化過程，主要包括邊界積分法和基于N-S方程與界面追蹤方法的多相流模擬方法。研究者們利用此類方法分析了液滴飛濺的臨界條件，指出液冠半徑與無量綱液冠演化時(shí)間成正比，其破碎的原因是液冠末端的Rayleigh-Taylor不穩(wěn)定性[7-10]。格子玻爾茲曼方法(lattice Boltzmann method，LBM)作為一種宏觀數(shù)值模擬方法，已經(jīng)發(fā)展成為研究復(fù)雜多相流現(xiàn)象的有效數(shù)值模擬方法，被用于多種液滴撞擊液膜的復(fù)雜多相流現(xiàn)象的模擬。與傳統(tǒng)宏觀多相流數(shù)值模擬方法相比，LBM優(yōu)勢(shì)[11-13]在于：①可以實(shí)現(xiàn)完全并行和分塊計(jì)算；②避免了傳統(tǒng)宏觀方程中的非線性對(duì)流項(xiàng)的處理，只需要求解一系列線性方程組；③壓力僅與流體狀態(tài)方程相關(guān)，不需要求解泊松方程來獲得壓力。其中LBM偽勢(shì)模型由于其簡(jiǎn)潔性和計(jì)算便利性被廣泛應(yīng)用，該模型不需要進(jìn)行界面跟蹤，通過粒子間作用的偽勢(shì)可以自動(dòng)形成兩相交界面[14]，已被成功用于液滴撞擊靜止液膜的研究[15-17]。目前基于液膜靜止過程的研究還缺乏對(duì)液冠演化過程的精確驗(yàn)證，無法保證模型模擬的準(zhǔn)確性[15-20]。由于液滴撞擊液膜過程中液膜往往無法保持靜止，因此研究液滴撞擊移動(dòng)液膜的過程具有非常重要的意義。

本文通過LBM偽勢(shì)模型模擬了液氣密度比為720條件下液滴撞擊移動(dòng)液膜的過程，分析了液膜運(yùn)動(dòng)速度等參數(shù)對(duì)液冠形態(tài)演化的影響，并將結(jié)果與現(xiàn)有半經(jīng)驗(yàn)公式進(jìn)行對(duì)比，以分析液冠半徑及液冠高度隨時(shí)間的變化規(guī)律，可為深入揭示復(fù)雜條件下液滴撞擊液膜過程相應(yīng)機(jī)理提供參考。

1 數(shù)值模型

本文主要采用多松馳時(shí)間(MRT)碰撞算子實(shí)現(xiàn)流體對(duì)流過程模擬，以獲得更好的數(shù)值穩(wěn)定性。采用MRT碰撞算子的相應(yīng)粒子分布方程[21]可以表示為

(1)

(2)

宏觀密度ρ和流場(chǎng)實(shí)際速度u可通過以下公式[11]求得

(3)

本文采用的外力格式為L(zhǎng)i等[21]提出的滿足熱力學(xué)一致性的外力格式

Fx,-Fx,Fy,-Fy,2(uxFx-uyFy),uxFx+uyFy)T

(4)

式中：ε為調(diào)節(jié)熱力學(xué)一致性的參數(shù)；ψ為粒子間的偽勢(shì)；Fx、Fy、ux、uy分別為作用在粒子上的合力F和流場(chǎng)實(shí)際速度u在x、y方向上的分量；Fm為粒子間的相互作用力。本文計(jì)算偽勢(shì)ψ時(shí)引入Carnahan-Starling(C-S)狀態(tài)方程[22]：

(5)

式中：a、b、R為相關(guān)參數(shù)，取值為a=0.25，b=4，R=1；Tc為臨界溫度；Pc為臨界壓力。

為獲得可單獨(dú)調(diào)節(jié)表面張力的模型，式(1)中源項(xiàng)C[21]為

(6)

式中Q可以通過下式求解[21]：

(7)

式中：κ為表面張力調(diào)節(jié)參數(shù)；G為粒子之間的相互作用強(qiáng)度；ωi為不同離散速度方向的權(quán)重系數(shù)。

2 模型驗(yàn)證

2.1 邊界及初始條件

液滴撞擊移動(dòng)液膜計(jì)算示意圖及液冠形狀演化相關(guān)參數(shù)如圖1所示，計(jì)算區(qū)域?yàn)? 001×301矩形區(qū)域，圖中r0為液滴初始半徑，H為液膜厚度，u0為液滴撞擊液膜的速度，uf為初始液膜運(yùn)動(dòng)速度。計(jì)算域上部邊界及左右邊界為非平衡外推邊界，底部邊界為無滑移邊界。結(jié)合液相密度ρl、液相黏滯系數(shù)υl、表面張力σ和時(shí)間t可以獲得液滴撞擊液膜時(shí)液冠演化形態(tài)的控制性無量綱參數(shù)：

(8)

式中：Re為雷諾數(shù)；We為韋伯?dāng)?shù);h*為液膜無量綱厚度;t*為無量綱時(shí)間;u*為初始無量綱液膜運(yùn)動(dòng)速度。液冠形狀演化過程中相應(yīng)參數(shù)如圖2所示，圖中rc為液冠半徑，Hu為液冠上游高度，Hd為液冠下游高度，對(duì)3個(gè)參數(shù)進(jìn)行無量綱化：

(9)

圖1 液滴撞擊液膜計(jì)算區(qū)域示意圖

圖2 液冠演化過程示意圖

液滴周邊和液膜附近密度場(chǎng)分別通過式(10)(11)初始化：

(10)

(11)

式中：(x0,y0)為液滴中心；ρg為初始汽相密度；y1為氣液交界面位置；w為初始界面寬度，考慮到計(jì)算穩(wěn)定性和物理真實(shí)性，設(shè)為5。

2.2 液滴撞擊靜止液膜驗(yàn)證

首先利用液滴撞擊靜止液膜的演化過程來驗(yàn)證模型的可靠性。液滴初始半徑r0=50，液膜厚度H=25。本文采用等溫模型，選取初始溫度T=0.5Tc，對(duì)應(yīng)液氣密度比ρl/ρg=720。液滴撞擊液膜速度設(shè)置為u0=0.125，初始無量綱液膜運(yùn)動(dòng)速度u*=0，松弛系數(shù)τυ=0.537 5，表面張力調(diào)節(jié)參數(shù)κ=0.2，對(duì)應(yīng)表面張力σ=0.008 3。相應(yīng)的Re和We分別為500和87.8。液滴撞擊液膜后液冠形態(tài)演化過程如圖3所示。

圖3 液滴撞擊靜止液膜不同時(shí)刻液冠形態(tài)演化過程

圖4 液滴半徑演化過程驗(yàn)證

圖5 液冠垂向延伸率演化過程驗(yàn)證

液滴半徑演化過程本文模型模擬結(jié)果與Gao等[6]半經(jīng)驗(yàn)公式計(jì)算結(jié)果對(duì)比如圖4所示，可見兩者結(jié)果一致。對(duì)比發(fā)現(xiàn)Gao等[6]模型試驗(yàn)結(jié)果在演化過程中，其半徑變化略小于公式計(jì)算值，其原因在于Gao等[6]在公式推導(dǎo)過程中假設(shè)在液冠形成之前液滴會(huì)與液膜完全融合，這一假設(shè)與試驗(yàn)現(xiàn)象[3]和數(shù)值模擬結(jié)果[10,14,18-20]存在差異。進(jìn)一步將本文模型模擬結(jié)果中的液冠垂向延伸率St與Gao等[6]半經(jīng)驗(yàn)公式計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果如圖5所示。在t*=0.500～1.200時(shí)，本文模型模擬結(jié)果與Gao等[6]半經(jīng)驗(yàn)公式計(jì)算結(jié)果吻合良好，但當(dāng)t*=1.200～1.500時(shí)，本文模型模擬結(jié)果突然增大，此時(shí)發(fā)生液冠不穩(wěn)定現(xiàn)象，液滴開始脫落，導(dǎo)致液冠垂向延伸率增大。隨著液滴脫落，在t*=1.500～2.000時(shí)，液冠垂向延伸率迅速下降，低于Gao等[6]半經(jīng)驗(yàn)公式計(jì)算結(jié)果。

3 模擬結(jié)果與分析

為分析不同Re、We和u*對(duì)液冠演化的影響，選取的相關(guān)參數(shù)如表1所示。模擬過程中保持液滴撞擊速度不變，通過調(diào)節(jié)τυ和κ獲得Re和We。

表1 液滴撞擊移動(dòng)液膜相關(guān)參數(shù)

3.1 Re的影響

不同Re條件下液滴撞擊液膜過程如圖6所示。隨著Re的增大，上、下游液冠厚度減小，同時(shí)液冠高度增大，液冠更容易發(fā)生破碎。受到初始無量綱液膜運(yùn)動(dòng)速度u*的影響，上、下游液冠演化過程呈現(xiàn)不對(duì)稱性。當(dāng)t*=1.250時(shí)，在Re=500和Re=1 000條件下上游液冠均已發(fā)生破碎；當(dāng)t*=2.000時(shí)，在Re=1 000條件下上游液冠已發(fā)生二次破碎，但在Re=200條件下上游液冠此時(shí)才發(fā)生破碎。不同Re條件下，整個(gè)演化過程中下游液冠并不發(fā)生破碎。

圖6 不同Re條件下液滴撞擊移動(dòng)液膜時(shí)液冠演化形態(tài)

圖7 不同Re條件下液滴撞擊移動(dòng)液膜時(shí)液冠參數(shù)演化過程

不同Re條件下，上、下游液冠高度和半徑演化過程如圖7所示。上、下游液冠高度均隨Re的增大而增大，且在同一時(shí)刻上游液冠高度大于下游液冠高度。此時(shí)液滴撞擊液膜后產(chǎn)生的環(huán)向速度與初始液膜運(yùn)動(dòng)速度相反，上游液冠處撞擊后液膜運(yùn)動(dòng)速度不連續(xù)性增大，促進(jìn)了上游液冠生長(zhǎng)和破碎過程。而下游液冠處撞擊后產(chǎn)生的環(huán)向速度與初始液膜運(yùn)動(dòng)速度方向相同，減小了下游液冠處撞擊后液膜運(yùn)動(dòng)速度不連續(xù)性，抑制了下游液冠的生長(zhǎng)。根據(jù)Gao等[6]的半經(jīng)驗(yàn)公式，在Re=200、500、1 000條件下對(duì)應(yīng)的能量損失系數(shù)分別為0.50、0.52和0.52，本文模型計(jì)算結(jié)果與Gao等[6]的半經(jīng)驗(yàn)公式計(jì)算結(jié)果十分吻合(圖7(c))。

3.2 We的影響

不同We條件下液滴撞擊液膜過程如圖8所示。上游液冠隨We的增大末端Rayleigh-Taylor不穩(wěn)定性增強(qiáng)，表面張力減小后液冠難以維持穩(wěn)定形態(tài)更易發(fā)生破碎。破碎液滴大小也隨We的增大而減小，在t*=1.250時(shí)，不同We條件下液冠均可以觀察到明顯的二次破碎。下游液冠穩(wěn)定性也隨著We的增大而減小，但均未發(fā)生破碎。當(dāng)We=1 165.5時(shí)下游液冠變形遠(yuǎn)大于We=87.8、139.4時(shí)下游液冠變形。

圖8 不同We條件下液滴撞擊移動(dòng)液膜時(shí)液冠演化形態(tài)

不同We條件下，上、下游液冠高度和液冠半徑演化過程如圖9所示。上游液冠高度不隨著We的增大而增大，而下游液冠高度隨著We的增大而增大。在同一時(shí)刻下游液冠高度大于上游液冠高度。上游液冠受到初始液膜運(yùn)動(dòng)速度的影響，Rayleigh-Taylor不穩(wěn)定性增強(qiáng)，末端極易發(fā)生破碎，導(dǎo)致上游液冠高度接近。但初始液膜運(yùn)動(dòng)速度抑制了下游液冠的Rayleigh-Taylor不穩(wěn)定性，致使下游液冠并不易發(fā)生液滴脫落。同時(shí)表面張力隨著We的增大而減小，液冠抗變形能力降低，大We時(shí)下游液冠變形更大，因此液冠高度更高。由圖9(c)可知，We的變化也不影響液冠半徑的演化。當(dāng)We=87.8、139.4、1 165.6時(shí)對(duì)應(yīng)的能量損失系數(shù)分別為0.52、0.51和0.49，We對(duì)液滴撞擊液膜過程中能量損失系數(shù)變化的影響可以忽略，本文模型計(jì)算結(jié)果與Gao等[6]的半經(jīng)驗(yàn)公式計(jì)算結(jié)果十分吻合。

圖9 不同We條件下液滴撞擊移動(dòng)液膜時(shí)液冠參數(shù)演化過程

3.3 初始無量綱液膜運(yùn)動(dòng)速度的影響

不同初始無量綱液膜運(yùn)動(dòng)速度條件下液滴撞擊液膜過程如圖10所示。上、下游液冠高度均隨初始無量綱液膜運(yùn)動(dòng)速度u*的增大而減小。當(dāng)t*=2.000時(shí)，在u*=0.8和u*=1.0條件下液冠末端均出現(xiàn)二次破碎。隨著u*增大，上游液冠處液滴撞擊產(chǎn)生的環(huán)向速度與初始液膜運(yùn)動(dòng)速度方向相反，上游液冠處撞擊后液膜運(yùn)動(dòng)速度不連續(xù)性增大，導(dǎo)致上游液冠末端Rayleigh-Taylor不穩(wěn)定性增強(qiáng)。而下游液冠處液滴撞擊產(chǎn)生的環(huán)向速度與初始液膜運(yùn)動(dòng)速度方向相同，下游液冠處撞擊后液膜運(yùn)動(dòng)速度不連續(xù)性減小，下游液冠Rayleigh-Taylor不穩(wěn)定性隨之減小。

圖10 不同初始無量綱液膜運(yùn)動(dòng)速度條件下液滴撞擊移動(dòng)液膜時(shí)液冠演化形態(tài)

不同初始無量綱液膜運(yùn)動(dòng)速度條件下上、下游液冠高度和液冠半徑演化過程如圖11所示。上、下游液冠高度均隨著初始無量綱液膜運(yùn)動(dòng)速度u*的增大而減小，上游液冠飛濺角受到初始液膜運(yùn)動(dòng)速度和液滴撞擊速度的共同作用，隨著u*增大，上游液冠飛濺角減小。液冠半徑隨u*增大而增大，當(dāng)u*=0.5、0.8、1.0時(shí)對(duì)應(yīng)的能量損失系數(shù)分別為0.52、0.58和0.62，隨著u*增大，撞擊過程中能量損失系數(shù)增大，由圖11(c)可知，本文模型計(jì)算結(jié)果與Gao等[6]的半經(jīng)驗(yàn)公式計(jì)算結(jié)果十分吻合。

圖11 不同初始無量綱液膜運(yùn)動(dòng)速度條件下液滴撞擊移動(dòng)液膜時(shí)液冠參數(shù)演化過程

4 結(jié) 論

a. 受到初始液膜運(yùn)動(dòng)速度的影響，上游液冠處撞擊后液膜運(yùn)動(dòng)速度不連續(xù)性增大， Rayleigh-Taylor不穩(wěn)定性加劇，更易發(fā)生破碎。而下游液冠處撞擊后液膜運(yùn)動(dòng)速度不連續(xù)性減小，下游液冠更穩(wěn)定。

b. 上、下游液冠高度均隨Re的增大而增大，但液冠半徑演化過程不隨Re的增大而改變。

c.We的變化影響到液冠末端穩(wěn)定性，上游液冠高度不隨We增大而改變，下游液冠隨We增大而發(fā)生劇烈變形，液冠高度隨We增大而增大。

d. 隨著初始無量綱液膜運(yùn)動(dòng)速度u*的增大，液滴撞擊過程中能量損失系數(shù)增大，上、下游液冠高度均減小，液冠半徑則隨u*的增大而增大。