液滴傾斜撞擊液膜液冠演化規(guī)律

郭平措，何小瀧，袁 浩

(1.青海省水利廳引黃濟(jì)寧工程建設(shè)管理局，西寧 810001； 2.四川大學(xué) 水力學(xué)及山區(qū)河流開發(fā)保護(hù)國家重點實驗室， 成都 610065； 3.重慶交通大學(xué) 重慶西南水運(yùn)工程科學(xué)研究所，重慶 400074)

1 研究進(jìn)展

液滴撞擊液膜現(xiàn)象通常伴隨著液滴撞擊干燥壁面，壁面液滴聚集形成液膜等現(xiàn)象，這一系列現(xiàn)象廣泛存在于工業(yè)過程和自然現(xiàn)象中，如泄洪霧化中二次激濺、農(nóng)藥噴灑等，因此液滴撞擊液膜過程中動力學(xué)機(jī)理的研究十分重要。

前人研究表明[1-3]，液滴撞擊液膜時液冠的演化過程被諸多參數(shù)所影響，包括氣液密度比、氣液運(yùn)動黏滯系數(shù)比、液膜厚度、撞擊速度、液相理化性質(zhì)等。當(dāng)前大部分研究主要通過實驗和數(shù)值模擬研究液滴垂直撞擊液膜時液冠的演化規(guī)律。Zhang[4]通過系統(tǒng)研究總結(jié)出液滴撞擊液膜時根據(jù)液冠的不同演化方式可以將碰撞過程歸納為融合、反彈、局部反彈、液冠、飛濺5種形態(tài)。Yarin等[5]通過實驗和邊界積分法研究了單個液滴撞擊液膜時液冠演化的規(guī)律。他指出液膜中的速度不連續(xù)性是液冠形成的主要原因，且液冠無量綱半徑與無量綱時間平方根之間存在線性關(guān)系。Gao等[6]通過大量實驗研究并結(jié)合理論分析給出了液滴撞擊靜止和移動液膜時液冠半徑演化的半經(jīng)驗公式。實驗研究可以通過高速攝像機(jī)準(zhǔn)確捕捉液冠的演化過程，但無法給出液冠演化過程中諸如流速場等內(nèi)部機(jī)理的演化過程。

與實驗研究相比，數(shù)值模擬可以有效獲得液滴撞擊液膜過程中速度場、渦量場等。通過數(shù)值模擬研究，Roisman等[7]和 Nikolopoulos等[8]指出液冠破碎的原因是液冠末端的Plateau-Rayleigh不穩(wěn)定現(xiàn)象。近30年來格子玻爾茲曼方法(LBM)多相流模型快速發(fā)展，現(xiàn)有多相流模型可以分為顏色模型、偽勢模型、自由能模型和相場模型4類[9]。Lee等[10]首先利用自由能模型模擬了大密度比條件下單個液滴撞擊液膜的過程。但該模型將滑移步分解為碰撞前滑移步和碰撞后滑移步，且需要求解各宏觀參數(shù)的方向偏導(dǎo)數(shù)，計算量較大。Mukherjee等[11]在Lee等[10]模型基礎(chǔ)上引入多松弛(MRT)碰撞算子，研究了密度比、黏滯系數(shù)比對液冠演化的影響。Cheng等[12]則進(jìn)一步將Mukherjee[11]二維模型轉(zhuǎn)化為三維模型，模擬了液滴撞擊移動液面的過程。

LBM偽勢模型由于計算簡潔性和高效性已經(jīng)被廣泛用于諸如沸騰[9]、空化[13]等復(fù)雜多相流現(xiàn)象的模擬中。相比于基于Navier-Stokes方程和界面捕捉方法的傳統(tǒng)數(shù)值模擬方法，該方法優(yōu)勢在于氣液交界面由粒子之間偽勢自動形成，不需要利用界面捕捉模型獲得；壓力與流體狀態(tài)方程相關(guān)，不需通過泊松方程求解壓力；僅需要求解一系列線性方程組，避免了對非線性對流項的處理[13]。Kharmiani等[14]和Yuan等[15]均利用LBM偽勢模型研究了單液滴垂直撞擊液膜時不同參數(shù)對液冠演化的影響。

當(dāng)前液滴撞擊液膜過程的研究主要集中于液滴垂直撞擊靜止或運(yùn)動液膜。因此本文通過大密度比LBM偽勢模型模擬了不同碰撞角條件下液冠演化過程，進(jìn)一步分析了雷諾數(shù)、韋伯?dāng)?shù)、液膜厚度、碰撞角度等對液冠演化的影響，為揭示復(fù)雜條件下液滴撞擊液膜過程打下基礎(chǔ)。

2 模型概述

液滴傾斜撞擊液膜計算示意圖及液冠形狀示意圖如圖1(a)所示，計算區(qū)域為1 001 lu×301 lu矩形區(qū)域。計算域左右邊界和頂部邊界均為非平衡外推邊界，底部為無滑移邊界條件[16]。圖中r0為液滴半徑；H為液膜厚度；U0為液滴撞擊液膜速度；θ為液滴撞擊速度與水平面之間夾角。結(jié)合液相密度ρl、液相運(yùn)動黏滯系數(shù)υl、表面張力σ和碰撞時間t，可獲得液滴撞擊液面時控制液冠形態(tài)變化的相應(yīng)參數(shù)，包括[15]

(1)

式中:Re為雷諾數(shù);We為韋伯?dāng)?shù);h*為液膜無量綱厚度;t*為無量綱時間。液冠形狀演化過程中相應(yīng)參數(shù)如圖1(b)所示，圖中S為液冠延伸長度，Hu為上游液冠高度;Hd為下游液冠高度。為便于進(jìn)行統(tǒng)一分析，本研究也對以上3個參數(shù)進(jìn)行無量綱化，即

(2)

圖1 液滴傾斜撞擊液膜過程計算及形狀演化示意圖Fig.1 Schematic diagram of the computation domain and the crown geometry

液滴周邊密度場通過式(3)進(jìn)行初始化：

(3)

式中：(x0,y0)為液滴中心；ρg為初始?xì)庀嗝芏龋粀為初始?xì)庖航唤缑鎸挾?。液膜附近流場初始化如?4)所示。

式中y1為氣液交界面位置。根據(jù)Li等[16]的研究，考慮到計算穩(wěn)定性和物理真實性，界面寬度w設(shè)置為lu。

3 LBM偽勢模型

為保證大密度比多相流現(xiàn)象模擬過程的數(shù)值穩(wěn)定性和物理真實性，本研究主要采用與MRT碰撞算子結(jié)合的偽勢模型。采用MRT碰撞算子粒子演化方程可以表示為[17]

(5)

其中ωi為相應(yīng)權(quán)重，本研究采用D2Q9格子模型，當(dāng)i=0時權(quán)重ωi=4/9，當(dāng)i=1～4時權(quán)重ωi=1/9，當(dāng)i=5～8時權(quán)重ωi=1/36。松弛矩陣Λ由不同松弛參數(shù)構(gòu)成[17]，即

(7)

(8)

宏觀密度ρ和流場實際速度u可通過對介觀粒子分布函數(shù)統(tǒng)計獲得[17],即

(9)

式中F為作用在粒子上的合力。為滿足模擬過程中的熱力學(xué)一致性，本研究采用Li等[18]提出的外力格式,即

(10)

式中：ε為調(diào)節(jié)熱力學(xué)一致性的參數(shù);ψ為粒子間的偽勢；Fm為粒子間的相互作用力，可以通過以下方程求解[18]

(11)

式中：G為粒子間的相互作用強(qiáng)度，wi為不同方向的權(quán)重，對于D2Q9模型其中i=0時權(quán)重系數(shù)wi=0，當(dāng)i=1～4時權(quán)重系數(shù)wi=1/3，當(dāng)i=5～8時權(quán)重系數(shù)[19]。Yuan等[20]將非理想流體狀態(tài)方程引入，獲得求解偽勢ψ的公式為

(12)

式中P為壓力。本研究采用Carnahan-Starling (C-S)狀態(tài)方程[21]

(13)

同時為使偽勢模型中的表面張力獨立于氣液密度比可進(jìn)行單獨調(diào)節(jié)，引入源項C[18]，即

(14)

式中Q可以通過式(15)求解[18]，即

(15)

式中κ為表面張力調(diào)節(jié)參數(shù)。

4 模型驗證

本研究首先通過液滴垂直撞擊靜止液膜過程來驗證模型的可靠性。本模擬中液滴半徑設(shè)置為r0=50 lu，液膜厚度設(shè)置為H=25 lu。模型采用等溫模型，初始溫度為T=0.5Tc，對應(yīng)液氣密度比為ρl/ρg=720。液滴初始碰撞速度設(shè)置為U0=0.125 lu·tu-1。松弛系數(shù)設(shè)置為τυ=0.537 5，表面張力調(diào)節(jié)參數(shù)設(shè)置為κ=0.2，對應(yīng)表面張力為σ=0.008 3。因此相應(yīng)的雷諾數(shù)和韋伯?dāng)?shù)分別為Re=500和We=87.8。

Yarin等[5]實驗研究表明，液滴撞擊靜止液膜過程中，無量綱液冠半徑與無量綱時間平方根成正比。前人研究中均為定性化驗證模型的可靠性[14-15,21-22]，本文則將LBM計算結(jié)果與理論分析結(jié)果進(jìn)行定量對比。Gao等[6]通過理論分析與大量實驗結(jié)果得到了準(zhǔn)確的液冠半徑r*與無量綱時間t*平方根之間的關(guān)系為

式中:r*為液冠中心到液冠射流底部無量綱半徑;β為與能量損失系數(shù)λs相關(guān)的參數(shù)[6]

(17)

其中能量損失系數(shù)λs與Re數(shù)、We數(shù)和無量綱液膜厚度h*之間相關(guān)[6]。

(18)

不同時刻液滴演化形態(tài)如圖2所示。液滴撞擊液面后液冠半徑與Gao等[6]提出的半經(jīng)驗公式對比如圖3所示。對比發(fā)現(xiàn)LBM計算結(jié)果與Gao等[6]數(shù)值模擬結(jié)果在液滴破碎前演變規(guī)律一致，而破碎后液冠半徑演化過程存在差異的原因是，Gao等[6]的模型未考慮到液冠破碎引起的動量變化。

圖2 不同時刻液冠演化形狀Fig.2 Crown geometry at different dimensionless instances

圖3 液冠半徑演化過程LBM數(shù)值模擬結(jié)果與理論 分析結(jié)果對比Fig.3 Comparison of the evolution process of crown radius between LBM numerical result and theoretical analysis

5 液滴傾斜撞擊液膜過程

為研究不同參數(shù)對液滴傾斜撞擊液膜時液冠演化過程的影響，本研究保持液滴初始碰撞速度U0=0.125 lu/tu不變，通過調(diào)節(jié)松弛系數(shù)τυ和表面張力調(diào)節(jié)系數(shù)κ獲得不同Re數(shù)與We數(shù)。不同參數(shù)如表1所示。

表1 液滴撞擊移動液膜相關(guān)參數(shù)Table 1 Parameters of a droplet obliquely impacts on liquid film

5.1 Re數(shù)影響

由于液滴傾斜碰撞液膜，導(dǎo)致上下游液冠演化出現(xiàn)不對稱性。圖4為無量綱時間t*=2.0時，不同Re數(shù)條件下液冠形態(tài)。Gao等[6]研究表明碰撞中液相損失能量隨著Re數(shù)增加而減小。當(dāng)Re=200時，上下游液冠均未發(fā)生破壞。但隨著Re數(shù)增大到500，碰撞中損失的能量減小，導(dǎo)致液膜中的流速不連續(xù)性增加，下游液冠發(fā)生破碎，伴隨著小液滴的形成，而上游液冠仍保持穩(wěn)定。

圖4 無量綱時間t*=2.0不同Re數(shù)條件下液冠演化 形態(tài)對比Fig.4 Crown geometry under varied Re number at dimensionless time t*=2.0

不同Re數(shù)條件下上下游液冠高度和液冠延伸長度演化如圖5所示。上游液冠高度隨著Re數(shù)的增加而增加，而下游液冠高度在液冠未發(fā)生破碎前也隨著Re數(shù)的增加而增加。隨著Re數(shù)增加，下游液冠更容易發(fā)生破碎。在下游液冠發(fā)生破碎后，Re=200和Re=500兩種雷諾數(shù)條件下下游液冠高度變化一致，而Re=1 000條件下同一時刻下游液冠高度更高。液冠延伸長度隨著時間變化如圖5(c)所示，由于大Re數(shù)條件下液滴撞擊液膜過程中消耗的能量更小，液冠未發(fā)生破碎前，液冠延伸長度隨著Re數(shù)的增加而增加。在下游液冠發(fā)生破碎后，3種Re數(shù)條件下液冠演化過程幾乎一致。

圖5 不同Re數(shù)條件下液冠演化過程Fig.5 Crown evolution process with varied Re number

5.2 We數(shù)影響

Kharmiani等[14]研究表明液冠演化過程中，隨著We數(shù)的增大，表面張力消耗的動能減小，導(dǎo)致液冠更不穩(wěn)定。圖6展示了不同We數(shù)條件下液冠在t*=2.0時刻的演化形態(tài)。隨著We數(shù)增加，液冠末端變得更薄和不穩(wěn)定。當(dāng)We=87.8時，上游液冠保持穩(wěn)定，而下游液冠發(fā)生一次破碎。隨著We數(shù)增加到139.4，上下游液冠均不能保持穩(wěn)定，發(fā)生破碎。而下游液冠破碎后生成的液滴尺寸小于上游液冠，破碎次數(shù)也隨著We數(shù)的增加而增加。

圖6 無量綱時間t*=2.0不同We數(shù)條件下液冠 演化形態(tài)對比Fig.6 Crown geometry under varied We number at dimensionless time t*=2.0

不同We數(shù)條件下上游液冠高度、下游液冠高度和液冠延伸長度演化如圖7所示。隨著We數(shù)增大，表面張力減小，液冠末端維持穩(wěn)定性能力降低，上游液冠高度在破碎前隨著We數(shù)的增加而增加。由于傾斜沖擊加強(qiáng)了下游液冠不穩(wěn)定性，因此下游液冠較上游液冠更早發(fā)生破碎，液冠破碎后高度仍隨著We數(shù)的增加而增加。圖7(c)對比了不同We數(shù)條件下液冠延伸長度隨時間的演化過程。Raman等[21]研究表明，We數(shù)對液冠延伸長度的影響可以忽略。本研究模擬結(jié)果表明，盡管隨著We數(shù)增大上下游液冠更容易發(fā)生破碎，但并不影響液冠延伸長度的演化，不同We數(shù)條件下液冠延伸長度一致。

圖7 不同We數(shù)條件下液冠演化過程Fig.7 Crown evolution process with varied We number

5.3 無量綱液膜厚度h*影響

Yuan等[15]研究表明，液膜厚度增加會導(dǎo)致液膜飛濺角增大和碰撞消耗能量增加。Gao等[6]更進(jìn)一步指出液膜厚度對碰撞中消耗的能量的影響遠(yuǎn)大于Re數(shù)和We數(shù)。圖8為在t*=2.0時刻液滴撞擊不同無量綱液膜厚度h*液冠形態(tài)。當(dāng)h*=0.1時，液膜厚度較小，碰撞中消耗能量較低，導(dǎo)致液冠末端更不穩(wěn)定，上下游液冠均發(fā)生破碎，伴隨著小液滴生成。同時由于液膜厚度較小，無法限制液冠射流的方向，液冠飛濺角較小。隨著液膜厚度增加到h*=0.25，此時上游液冠均保持穩(wěn)定，而下游液冠發(fā)生破碎。隨著液膜厚度增加到h*=0.5，上下游液冠均保持穩(wěn)定，而液冠飛濺角隨著液膜厚度的增加而增加。

圖8 無量綱時間t*=2.0不同無量綱液膜厚度h* 條件下液冠形態(tài)Fig.8 Crown geometry at t*=2.0 with varied dimensionless liquid film thickness h*

圖9 不同無量綱液膜厚度h*條件下液冠演化過程Fig.9 Crown evolution process with varied dimensionless liquid film thickness h*

圖9展示了不同無量綱液膜厚度h*條件下，上游液冠高度、下游液冠高度和液冠延伸長度隨時間變化過程。在液冠破碎前，無量綱液膜厚度由h*=0.1增加到h*=0.25，受到液冠飛濺角的影響，上游液冠高度隨液膜厚度增加而增加。但無量綱液膜厚度由h*=0.25增加到h*=0.5，隨著液膜厚度增加，碰撞中損失能量增加，上游液冠高度反而降低。下游液冠演化過程與上游液冠相似。Gao等[6]研究表明，液滴撞擊液膜過程中消耗能量越多，液冠半徑演化速率越慢。不同無量綱液膜厚度h*條件下液冠延伸長度隨時間變化規(guī)律如圖9(c)所示，盡管h*=0.1時，上下游液冠發(fā)生多次破碎，但其液冠半徑演化速率大于其余兩種液膜厚度。

5.4 碰撞角θ影響

為進(jìn)一步研究碰撞角對液冠演化的影響，本文保持液滴速度不變，改變液滴初始速度與水平方向的夾角。t*=2.0時刻不同碰撞角θ條件下液冠形態(tài)如圖10所示。在θ=45°時，下游液冠并未發(fā)生破碎，但隨著碰撞角增加到θ=60°，液冠發(fā)生破碎。上游液冠根部厚度則隨著碰撞角θ增大而減小，上下游液冠高度均增大。

圖10 無量綱時間t*=2.0不同撞擊角θ條件下液冠 演化形態(tài)對比Fig.10 Crown geometry under varied collision angle θ at dimensionless time t*=2.0

圖11展示了不同碰撞角θ條件下，上游液冠高度、下游液冠高度和液冠延伸長度隨時間變化過程。碰撞時液滴垂向速度隨著θ的增大而增大，導(dǎo)致碰撞過程中能量損失減小。上游液冠高度、下游液冠高度及液冠延伸長度均隨著碰撞角θ的增大而增大，且下游液冠高度隨著碰撞角θ的增大更早發(fā)生破碎。

圖11 不同撞擊角θ條件下液冠演化過程Fig.11 Crown evolution process with varied impact angle θ

6 結(jié)論與展望

本研究利用可調(diào)節(jié)表面張力的大密度比LBM偽勢模型，研究了液滴傾斜撞擊液膜的過程，并分析了不同Re數(shù)、We數(shù)、無量綱液膜厚度h*和碰撞角θ對液冠演化過程的影響，得到以下結(jié)論：

(1)碰撞時能量損失隨著Re數(shù)增大而減小，進(jìn)而導(dǎo)致上下游液冠高度和液冠延伸長度隨著Re數(shù)的增大而增大。

(2)隨著We數(shù)增大，上下游液冠更易發(fā)生變形和破壞，但不同We數(shù)條件下液冠延伸長度接近。

(3)液膜厚度增加會導(dǎo)致碰撞時消耗的能量更多，并增大液冠飛濺角。無量綱液膜厚度存在一個閾值h*=0.25，使未破碎的上下游液冠高度達(dá)到最大值，液冠延伸長度隨著液膜厚度的增加而減小。

(4)液滴傾斜撞擊液面會導(dǎo)致液冠的不對稱發(fā)育，上下游液冠高度和液冠延伸長度均隨著碰撞角θ的增大而增大。

本研究初步探討了二維液滴傾斜撞擊液膜的特性，考慮到液滴傾斜液膜會造成液冠的不對稱性和復(fù)雜性，二維模型難以模擬出液冠中流速的演化過程，無法深入分析獲得液冠上下游發(fā)生不對稱演化和破碎的機(jī)理。因此在后續(xù)研究中，本文所使用LBM偽勢模型將被拓展到D3Q19格子模型上，并對復(fù)雜壁面條件下液滴撞擊液膜后液冠的演化過程進(jìn)行深入研究。