微小水滴撞擊深水液池空腔運動的數值模擬及機理研究?

裴傳康 魏炳乾

(西安理工大學,省部共建西北旱區生態水利國家重點實驗室,西安 710048)

(2018年7月28日收到;2018年9月20日收到修改稿)

1 引 言

液滴下落沖擊不同介質的運動過程在科技應用及自然現象中廣泛存在.液滴沖擊液體表面的研究在犯罪取證、開發不可浸潤表面或完全可浸潤表面、高精度活化或表面污染物轉移等方面有諸多應用[1],液-液界面的部分聚結過程也是許多復雜物理現象的組成部分,與地球物理學與土力學等有諸多關聯[2,3];對于海洋、湖泊等大面積水面,其廣泛的自然摻氣現象主要取決于水滴撞擊所引發的氣泡夾帶[4];在水利工程中,高水頭泄水建筑物的霧化、摻氣、消能等問題也與此息息相關.因此,研究此類基本運動對于理解自然界氣液流動的界面變形、改善液滴運動在工程中的應用具有重要的意義與作用.

1963年,Worthington[5]通過實驗首次系統地描述了液滴和固體小球撞擊液池的過程及其運動規律.隨后的大量研究表明,當液滴以較低的撞擊速度撞擊液面時,液滴與液面發生完全聚結現象,并在液面下方生成渦環[6,7],在這種情況下,自由液面并不會出現飛濺射流,而是呈現出平坦的狀態.隨著液滴撞擊速度的增加,自由液面開始飛濺,液面出現較大變形,產生一個空腔,并在四周形成皇冠狀的射流[8],空腔塌陷后,液柱從撞擊中心升起,速度較大時,由于Plateau-Rayleigh不穩定性的影響,會在液柱上方分離出一個或多個小液滴[9].一般情況下,液滴沖擊液池的運動主要受重力、慣性力與表面張力的影響,因此對于弗勞德數(Fr)和韋伯數(We)十分敏感[10].

傳統上對液滴運動的探索多是以高速攝像機為主要工具進行的實驗研究,但由于液滴運動的復雜性,實驗攝影方式較難捕捉細致的自由面邊界,對諸如壓力和速度等物理量的測定也十分不便.隨著Hirt與Nicolas在1981年[11]提出利用流體體積法(VOF)分離兩種流體構成的尖銳界面后,數值模擬在捕捉界面變形方面精度不斷提升[12,13],逐漸成為了研究液滴撞擊運動的有力工具.Yue等[14]對液滴撞擊液面運動進行了基于相場法的數值模擬,并根據最大Oh數得到了部分融合發生的臨界標準.Ray等[15]利用CLSVOF法對液滴撞擊液面進行了更為具體的研究,結果表明慣性力和表面張力在很大程度上決定了液面的運動過程.Castillo-Orozco等[16]基于VOF法和表面張力模型(CSF),討論了流體的物理化學性質和液滴沖擊速度變化對二次液滴的斷裂以及冠狀射流的影響.戴劍鋒等[17]應用CLSVOF法研究了液滴撞擊傾斜表面液膜后液膜的形態演化及飛濺過程,研究表明,空氣卷吸氣泡數量隨著沖擊角的增大而減小.黃虎等[18]則建立了液滴撞擊平面液膜的數學模型,采用格子Boltzmann方法探討了相對液膜厚度和表面張力等物理參數對界面運動過程的影響.

以往對液滴撞擊運動的研究多針對毫米級直徑的液滴[19],驅動其空腔形成的要素主要為慣性力及重力[10],但在較小液滴尺度下,由于液滴質量及體積的減小,毛細波運動在驅動液-液界面變形及空腔形成過程中起到越發重要的作用.現有研究較少關注直徑在微米層級的液滴運動及空腔形成過程以及決定微小液滴空腔運動的主要驅動力,因此,本文主要關注高速運動的微小水滴,運用基于四叉樹結構自適應網格、VOF方法以及變密度不可壓縮Navier-Stocks方程的開源多相流程序Gerris[20?22],研究一定物理參數范圍下微小水滴撞擊深水液池的液面變形及毛細波運動過程,并嘗試揭示水滴撞擊產生的空腔運動及形成機理.

2 計算區域及數值方法

2.1 控制方程

液滴撞擊液池的運動過程可以使用帶有表面張力項的變密度、不可壓縮Navier-Stocks方程來描述,具體控制方程如下[20,21]:

式中,ρ=ρ(x,t)為流體密度,u=(u,v,w)為流體速度,p為壓力,μ=μ(x,t)是流體的動力黏度,變形張量D定義為Dij=(?iuj+?jui)/2,σ為表面張力系數,κ為界面曲率,狄拉克分布函數δs表示表面張力僅作用于兩相界面處,n為兩相界面的法向量.

Gerris采用經典的VOF方法追蹤相界面,對于兩相流動,引入計算網格中第一種流體的體積分數c(x,t),并定義混合流體的密度和黏度為:

式中,ρ1,ρ2,μ1,μ2分別是第一種流體和第二種流體的密度以及黏度;函數?c由體積分數c平滑處理后得出,以便提高計算的穩定性.

密度對流方程可由等效的體積分數對流方程替換

2.2 數值方法

本文數值模擬采用基于Linux的開源軟件Gerris進行,該軟件使用基于四叉樹(二維)/八叉樹(三維)的自適應空間離散方法,使用分步投影方法求解變密度不可壓縮的Navier-Stocks方程,使用VOF方法跟蹤相界面.高度函數和界面附近的自適應網格細化可以精確表示表面張力作用,對流項使用Godunov格式求解,并行計算采用MPI庫進行.

如圖1所示,水滴撞擊深水液池的數值模擬在軸對稱坐標系中進行,Y軸為計算區域的對稱軸,D為初始水滴直徑,R=D/2,正方形計算區域的長度H=20D,水滴距離液面的距離H1=0.1D,液池深度H2=12D以消除底部對液滴撞擊運動的影響,水滴在重力g和撞擊速度Vi的作用下撞擊液池.采用雷諾數、韋伯數和弗勞德數來描述液滴撞擊的運動特征,三者分別表征液體慣性力與黏滯力間的關系、液體慣性力與表面張力間的關系以及液體慣性力與重力間的關系.三個無量綱參數的表達式如(7)式所示,主要物理參數如表1所列.

表1 主要物理參數Table 1.Main physical parameters.

圖1 計算區域簡圖Fig.1.Schematic diagram of the computational domain.

2.3 自適應網格技術

采用數值方法對高速液滴的撞擊運動進行準確的模擬極具挑戰性,因為運動產生的微小界面變形、復雜的幾何形狀以及特征尺度的巨大差異需要足夠的網格分辨率來捕捉,從而大幅地增加了計算量與計算時間.目前針對該問題的一個有效解決方法是采用自適應網格(adaptive mesh refinement,AMR)技術.根據流動特征對網格進行局部細化或粗化使得AMR技術可以將計算效率集中在最需要的區域,從而以最小的計算成本獲取精確的結果.

本文采用Gerris進行數值模擬,Gerris使用有限體積法(FVM)來求解控制方程,并根據四叉樹網格自適應規則和條件將計算域離散為不同等級的計算網格.水滴撞擊深水液池數值模擬的關鍵位置在于液-液界面的交接處以及相界面附近,本文依此設計如下網格自適應規則,每一步更新一次計算網格,其中最大網格加密層數為11層,即在一個計算區域(box,Lbox=10)內的最大網格數量為211.圖2為計算區域初始狀態自適應網格的空間離散示意.

1)計算初始加密水滴與液池接觸區域,即水滴與液池相界面處正負0.15內的網格至11層.

2)自動加密相界面附近體積分數在0—1之間、梯度變化劇烈區域的網格,最大加密到11層,最小加密到6層,以最小化界面重建產生的誤差.

3)自動加密渦量變化區域的網格,根據其變化劇烈程度最大加密到11層,最小加密到4層.

4)根據U,V速度分量的變化自動加密網格,最大加密到11層.

5)限制2)—4)條規則最小加密層數的加密區間為:Y向液滴中心上方2R至水面下方4R;x向對稱軸左右4R內的矩形內,以提高計算效率.

圖2 計算初始狀態的空間離散Fig.2.The spatial discretization of the initial simulation.

2.4 模型驗證與率定

為了保證數值模擬結果的準確性,本文選擇Morton的實驗數據[23]對數值模型進行驗證,實驗使用直徑為2.9 mm的液滴撞擊液池,弗勞德數及韋伯數分別為Fr=220,We=248.如圖3所示,照片為高速攝影機攝得的實驗過程,白色線條表示相同時間節點下的數值模擬結果,t為物理時間乘以Vi/D后的無量綱時間.液滴下落后沖擊液池并產生了一個空腔,腔體在t=7.9時達到最大化.空腔塌陷后毛細波向中心處傳遞,并坍縮形成中心射流,使其高度不斷增大,在射流頂端斷裂生成二次液滴.由于實驗環境的復雜性,模擬條件與實驗條件無法完全一致,且本文采用軸對稱模型假定進行模擬,無法捕捉非對稱運動,因此模擬值與實驗值存在一定差異.但數值模擬在界面變形、空腔的形成與成長、毛細波在空腔底部的傳播等方面與實驗值取得了良好的一致性,且在空腔形成過程中給出了較實驗更加詳盡的毛細波運動細節,中心射流最大高度以及空腔最大深度的誤差分別為1.7%,2.6%,表明數值模擬能夠較好地描述液滴撞擊液池的運動.

圖3 數值模擬與實驗結果對比Fig.3.Comparison of the numerical and experimental results.

3 計算結果與分析

3.1 基本撞擊過程

水滴撞擊深水液池后相界面的變化可以直觀反映出水滴撞擊的運動特性、水滴與池水的混摻過程以及撞擊所夾帶的氣泡大小.因此,使用經過驗證的數值方法研究直徑為290μm的水滴以2.5—6.5 m/s的速度撞擊深水液池的運動過程.

圖4 不同工況下自由液面隨時間的運動過程 (a)Fr=117.2,Re=725,We=25.2,Vi=2.5 m/s;(b)Fr=229.7,Re=1015,We=49.3,Vi=3.5 m/s;(c)Fr=379.7,Re=1305,We=81.6,Vi=4.5 m/s;(d)Fr=567.1,Re=1595,We=121.8,Vi=5.5 m/s;(e)Fr=792.1,Re=1885,We=170.2,Vi=6.5 m/sFig.4.Free surface profiles simulated at selected times:(a)Fr=117.2,Re=725,We=25.2,Vi=2.5 m/s;(b)Fr=229.7,Re=1015,We=49.3,Vi=3.5 m/s;(c)Fr=379.7,Re=1305,We=81.6,Vi=4.5 m/s;(d)Fr=567.1,Re=1595,We=121.8,Vi=5.5 m/s;(e)Fr=792.1,Re=1885,We=170.2,Vi=6.5 m/s.

水滴以五種不同撞擊速度進入水池后自由表面的模擬輪廓以及水滴分布如圖4所示.藍色為氣相,紅色為水池部分液相,天藍色為水滴部分液相,從工況(a)—(e),水滴的撞擊速度不斷增大.在圖4(a)中,Fr=117.2,We=25.2,Re=725,水滴以較低速度撞擊深水液池,發生完全聚結現象,并在水池中形成底部夾帶一個氣泡的空腔,水滴入水后基本附著在自由液面之下.隨著時間增加,空腔開始向內坍縮,同時毛細波向空腔底部傳遞,使得水滴部分的水體向中心聚合并產生兩個對稱的渦,最終在空腔塌陷后對池水產生穿刺效應,池水在慣性力的作用下逐漸平復.在圖4(b)—(d)中,撞擊產生的空腔隨著水滴撞擊時動能的增加越來越大,空腔底部夾帶的氣泡不斷變小,隨著撞擊速度的增大,在空腔塌陷后,一個短而粗的射流在撞擊中心產生,且射流最大高度不斷增加.由于毛細波的傳遞速度加快,水滴入水后附著自由液面的面積也逐漸變大,產生的渦距離撞擊中心越來越遠.在圖4(e)中,Fr=792.1,Re=1885,We=170.2,水滴撞擊后首先產生一個U形空腔,最大深度再次增加,隨后由于毛細波的快速傳播推動空腔底部部分側壁收縮閉合,截留形成一個較大的氣泡,并在閉合處產生極細的射流,射流由于豎向剪應力較大,引起Plateau-Rayleigh不穩定性,在尖端斷裂生成多個二次滴,自由液面變化更為劇烈.

上述水滴撞擊水池的過程可以歸納為三個基本階段,第一階段為空腔的形成與膨脹;第二階段為毛細波傳播導致的空腔收縮;第三階段為空腔的回復.液滴撞擊深水液池時的運動分為完全聚結與部分聚結現象[24],在水滴與深水液池完全聚結時,水滴入水時由于水滴底部發生凹陷變形產生的氣泡夾帶隨著撞擊速度的增加而減小,而在部分聚結發生時氣泡的夾帶與截留行為則更為復雜.

3.2 空腔運動的基本特性

為了深入探究微米級水滴撞擊深水液池后空腔運動的基本特性,對本文五個模擬工況進行定量研究.描述空腔幾何特征的基本參數如圖5所示,以初始靜水狀態下深水液池的液面高度為基準線,hw為基準線至波浪頂端的波浪隆起高度,h為基準線至空腔最低點的空腔深度,r為空腔基準線上軸對稱點處至空腔側壁的空腔寬度.下文涉及的所有長度單位均為實際長度除以初始水滴直徑D后的無量綱長度,時間t為實際時間乘以Vi/D后的無量綱時間.

圖6為不同工況下空腔深度h隨時間的變化過程.由圖6可以看出,液滴撞擊深水液池后,由于動能、空腔側壁隆起部分重力勢能以及表面張力能的驅動,空腔深度先以較快速度增長,其后增勢逐漸放緩,在到達最大空腔深度后快速回彈,至接近原自由液面后回彈速度逐漸放緩,近乎趨于直線.弗勞德數的增加對空腔深度變化影響顯著,在工況(a)中,弗勞德數僅為117.2,空腔能量耗散時程較短,回復較為迅速,隨著弗勞德數的增大,撞擊動能增大,回復時程也逐漸增加.值得注意的是,在工況(e)中,由于空腔形狀由U形轉變為近似梯形,空腔在t=1.5—2.6時先緩慢上升,隨后毛細波向空腔底部傳遞,促使其底部變為圓柱形,且深度繼續增加,最后空腔底部側壁坍塌形成射流,因此在t=2.6后空腔回復速度較其他工況更為迅速.將數值結果運用最小二乘法擬合,得到的擬合曲線表達式如(8)式所示,該式能夠在忽略毛細波運動的前提下,在空腔深度為h=D至達到最大深度的范圍內較好地描述了空腔深度隨時間的成長關系,Liow等[8]與Ray等[24]對液滴撞擊運動中時間與空腔深度的關系進行了建模,得出與本文相似的結論:

圖5 空腔的幾何特征示意Fig.5.Geometric characteristics of the cavity.

圖6 不同工況下空腔深度隨時間的變化Fig.6.Cavity depth with time under different conditions.

圖7及圖8分別為不同工況下波浪隆起高度hw和空腔寬度r隨時間的變化過程.由圖7易得,水滴撞擊水池后波浪高度的變化過程也經歷了從快速增長到逐漸回落,最終在慣性力的作用下自主運動,趨于穩定的過程.隨著弗勞德數的增大,更大的垂向速度在創造更深空腔的同時激發了更高的最大波浪隆起高度.從圖8可以看出,空腔寬度的增長與波浪的運動息息相關,在撞擊產生的動能與波浪自身重力勢能基本耗散后,空腔運動主要由毛細波及其干擾驅動,最后在慣性下線性增長.毛細波現象在圖8工況(e)的t=0.8時非常顯著,其周期在圖6的t=1.5—2.5中也有體現.

圖7 不同工況下波浪隆起高度隨時間的變化Fig.7.Wave height with time under different conditions.

圖8 不同工況下空腔寬度隨時間的變化Fig.8.Cavity width with time under different conditions.

3.3 空腔形成以及毛細波傳播機理

選取水滴完全聚結的工況(d)與水滴部分聚結的工況(e)對空腔形成以及毛細波傳播的機理進行深入研究.圖9為不同時間節點下工況(d)及工況(e)水滴撞擊深水液池空腔運動的等值線圖,其中左側為渦量場等值線圖,右側為壓力場等值線圖,黑色線條表示相界面.由圖9(a)可得,水滴撞擊液池后,空腔形狀由U形向V形轉變,最終腔底首先上升,形成沒有氣泡截留的射流.

Berberovic等[25]的研究表明,毛細波的傳播路徑與低壓區的形成息息相關.在圖9(a)壓力場等值線圖中可以看出,低壓區首先在空腔側壁與底部的交界處產生,隨后向空腔底部傳播,并形成一個相對尖銳的點,在空腔底部上升后,毛細波開始向下方傳遞,并逐漸趨于平緩.

在發生部分聚結現象的圖9(b)中,低壓區首先在波浪底部與側壁上交界處產生,并已經在相界面上形成了一個尖銳的折點,此時空腔形狀為半球形,隨后毛細波向空腔底部傳播,在底部中心空腔轉變為圓柱狀,低壓區附著在圓柱下方與底面交點,毛細波在空腔坍塌前并沒有到達空腔底部中心,而后空腔側壁坍塌形成氣泡截留,并在中心產生了快速射流液柱.

渦量定義為流體速度矢量的旋度,描述液體流動的剪切運動.在工況(d)中,如圖9(a)渦量場所示,流體渦量一直跟隨毛細波位置,當空腔坍塌產生慢速射流時,渦量在靠近液面區域以及空腔底部靠近對稱軸的區域各產生一個較大的渦環,而在工況(e)中,渦環的生成被抑制,流體僅在t=1.906時表現出一個靠近空腔底部的小渦環,而后渦環迅速消失.由渦量場與壓力場對比可得,渦量較大的區域并不總是處于低壓區內,撞擊運動初始自由液面的壓力差可能是渦量產生的誘因,但其后低壓區的運動與渦量之間并無較強的相關性.

4 結 論

本文采用基于自適應網格和VOF方法的開源程序Gerris對微小水滴撞擊深水液池后的流動過程以及空腔生長進行了數值模擬,研究了不同Fr數對撞擊后空腔毛細波運動以及氣泡截留的影響,主要得到以下結論.

1)在恰當的自適應條件下,Gerris程序能夠在節約計算資源的同時較為準確地預測水滴撞擊深水液池的運動,數值模擬所得的界面變形、空腔成長、毛細波的傳播以及中心射流過程與實驗結果符合良好.

圖9 不同時間節點下水滴撞擊深水液池空腔運動的渦量場和壓力場等值線圖 (a)Fr=567.1,Re=1595,We=121.8;(b)Fr=792.1,Re=1885,We=170.2Fig.9.Vorticity and pressure contours at selected time:(a)Fr=567.1,Re=1595,We=121.8;(b)Fr=792.1,Re=1885,We=170.2.

2)液滴下落撞擊深水液池后,液面擴張形成一個空腔,其后隨著毛細波運動逐漸回縮.液滴完全聚合時,空腔形狀往往由U形向V形轉變;在液滴部分聚合生成細長中心射流并產生氣泡截留時,空腔初始形狀則近似一個半球形,其后在底部變形為圓柱形.

3)液滴撞擊深水液池后,空腔深度先以較快速度增長,在到達最大空腔深度后快速回彈,至接近原自由液面后速度逐漸放緩.在忽略毛細波作用、空腔深度為h=D至h=hmax范圍內的前提下,空腔深度隨時間的成長關系可由(Vit)/D=0.15(h/D)5/2來描述,但最終空腔底部的形成是由毛細波運動決定的;空腔寬度的增長主要由毛細波運動及其干擾驅動,最后在慣性力作用下線性增長.

4)毛細波運動可由壓力場中低壓區的位置示蹤,在撞擊速度較低,液滴完全聚結時(Fr=567.1,Re=1595,We=121.8),低壓區首先在空腔側壁與底部交界處產生,隨后向空腔底部傳播,在靠近液面以及空腔底部靠近中心區域各產生一個較大的渦環;在發生部分聚結現象,產生細長射流時(Fr=792.1,Re=1885,We=170.2),渦環的生成被抑制,低壓區首先在波浪底部與側壁上交界處產生,空腔底部變為圓柱狀后,毛細波在空腔坍塌前并沒有到達空腔底部中心,導致空腔側壁首先坍塌形成氣泡截留.

感謝西安理工大學水力學研究所的李林博、楊泓、薛博升和楊琰青在論文完成過程中的幫助與討論.