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数值模拟与实验结果对比