2計(jì)算區(qū)域及數(shù)值方法

2.1控制方程

液滴撞擊液池的運(yùn)動(dòng)過(guò)程可以使用帶有表面張力項(xiàng)的變密度、不可壓縮Navier-Stocks方程來(lái)描述，具體控制方程如下：

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

Gerris采用經(jīng)典的VOF方法追蹤相界面，對(duì)于兩相流動(dòng)，引入計(jì)算網(wǎng)格中第一種流體的體積分?jǐn)?shù)c（x,t），并定義混合流體的密度和黏度為：

式中，ρ1,ρ2,μ1,μ2分別是第一種流體和第二種流體的密度以及黏度；函數(shù)?c由體積分?jǐn)?shù)c平滑處理后得出，以便提高計(jì)算的穩(wěn)定性。

密度對(duì)流方程可由等效的體積分?jǐn)?shù)對(duì)流方程替換

2.2數(shù)值方法

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

如圖1所示，水滴撞擊深水液池的數(shù)值模擬在軸對(duì)稱坐標(biāo)系中進(jìn)行，Y軸為計(jì)算區(qū)域的對(duì)稱軸，D為初始水滴直徑，R=D/2,正方形計(jì)算區(qū)域的長(zhǎng)度H=20D,水滴距離液面的距離H1=0.1D,液池深度H2=12D以消除底部對(duì)液滴撞擊運(yùn)動(dòng)的影響，水滴在重力g和撞擊速度Vi的作用下撞擊液池。采用雷諾數(shù)、韋伯?dāng)?shù)和弗勞德數(shù)來(lái)描述液滴撞擊的運(yùn)動(dòng)特征，三者分別表征液體慣性力與黏滯力間的關(guān)系、液體慣性力與表面張力間的關(guān)系以及液體慣性力與重力間的關(guān)系。三個(gè)無(wú)量綱參數(shù)的表達(dá)式如（7）式所示，主要物理參數(shù)如表1所列。

表1主要物理參數(shù)

圖1計(jì)算區(qū)域簡(jiǎn)圖

2.3自適應(yīng)網(wǎng)格技術(shù)

采用數(shù)值方法對(duì)高速液滴的撞擊運(yùn)動(dòng)進(jìn)行準(zhǔn)確的模擬極具挑戰(zhàn)性，因?yàn)檫\(yùn)動(dòng)產(chǎn)生的微小界面變形、復(fù)雜的幾何形狀以及特征尺度的巨大差異需要足夠的網(wǎng)格分辨率來(lái)捕捉，從而大幅地增加了計(jì)算量與計(jì)算時(shí)間。目前針對(duì)該問(wèn)題的一個(gè)有效解決方法是采用自適應(yīng)網(wǎng)格（adaptive mesh refinement,AMR）技術(shù)。根據(jù)流動(dòng)特征對(duì)網(wǎng)格進(jìn)行局部細(xì)化或粗化使得AMR技術(shù)可以將計(jì)算效率集中在最需要的區(qū)域，從而以最小的計(jì)算成本獲取精確的結(jié)果。

本文采用Gerris進(jìn)行數(shù)值模擬，Gerris使用有限體積法（FVM）來(lái)求解控制方程，并根據(jù)四叉樹網(wǎng)格自適應(yīng)規(guī)則和條件將計(jì)算域離散為不同等級(jí)的計(jì)算網(wǎng)格。水滴撞擊深水液池?cái)?shù)值模擬的關(guān)鍵位置在于液~液界面的交接處以及相界面附近，本文依此設(shè)計(jì)如下網(wǎng)格自適應(yīng)規(guī)則，每一步更新一次計(jì)算網(wǎng)格，其中最大網(wǎng)格加密層數(shù)為11層，即在一個(gè)計(jì)算區(qū)域（box,Lbox=10）內(nèi)的最大網(wǎng)格數(shù)量為211.圖2為計(jì)算區(qū)域初始狀態(tài)自適應(yīng)網(wǎng)格的空間離散示意。

1）計(jì)算初始加密水滴與液池接觸區(qū)域，即水滴與液池相界面處正負(fù)0.15內(nèi)的網(wǎng)格至11層。

2）自動(dòng)加密相界面附近體積分?jǐn)?shù)在0—1之間、梯度變化劇烈區(qū)域的網(wǎng)格，最大加密到11層，最小加密到6層，以最小化界面重建產(chǎn)生的誤差。

3）自動(dòng)加密渦量變化區(qū)域的網(wǎng)格，根據(jù)其變化劇烈程度最大加密到11層，最小加密到4層。

4）根據(jù)U,V速度分量的變化自動(dòng)加密網(wǎng)格，最大加密到11層。

5）限制2）—4）條規(guī)則最小加密層數(shù)的加密區(qū)間為：Y向液滴中心上方2R至水面下方4R;x向?qū)ΨQ軸左右4R內(nèi)的矩形內(nèi)，以提高計(jì)算效率。

圖2計(jì)算初始狀態(tài)的空間離散

2.4模型驗(yàn)證與率定

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

圖3數(shù)值模擬與實(shí)驗(yàn)結(jié)果對(duì)比