2計算區域及數值方法

2.1控製方程

液滴撞擊液池的運動過程可以使用帶有表麵張力項的變密度、不可壓縮Navier-Stocks方程來描述，具體控製方程如下：

式中，ρ=ρ（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主要物理參數

圖1計算區域簡圖

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計算初始狀態的空間離散

2.4模型驗證與率定

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

圖3數值模擬與實驗結果對比