基于VOF模型的结构物出水过程数值模拟

邹 星，李海涛，宗 智

(大连理工大学船舶工程学院，辽宁 大连 116024)

结构物出水在许多工程领域都具有非常重要的意义，一直是流体力学工作者研究的重要内容。当物体从水中上升，穿过自由液面并进入空气时，周围介质的密度会发生一个很大的阶跃，从水的密度998.2 kg/m3突然变为1.225 kg/m3，这种突变类似于气体动力学中的强间断问题，就像穿过了一层很强的激波一样[1]，在水面附近会发生一系列复杂的物理现象。

出水问题的理论研究可追溯到 KARMAN[2]的工作，他在1929年提出了计算物体撞击载荷的理论模型，这一思想至今不失其光彩。HAVELOCK[3]用浸没于水中的运动的偶极子来计算水平放置的圆柱体在自由表面下的水平运动，给出了瞬时的受力结果。LIJU等[4]利用边界元方法(BEM)模拟了一个对称刚体在垂直穿过自由面时的液面变化，并利用实验进行了验证。

国内学者在结构物出水方面也做了很多工作。鲁传敬[5]使用势流模型得到细长体出入水问题的解析解。蔡军辉等[6]通过解N-S方程，采用MAC法进行平头圆柱体垂直出水黏性流动的数值模拟。在实验方面，张军等[7]采用 TRPIV技术对钝头回转体垂直出水进行测量，揭示了出水过程中流动结构及其演变。

VOF模型通过求解单独的动量方程和处理穿过区域的每一流体的体积分数来模拟两种或3种不能混合的流体[8]。在VOF模型中，两种流体没有互相穿插，对增加到模型里的每一附加相，就引进一个变量，即计算单元里的相的体积分数[9]。在每个控制体积内，所有相的体积分数为1。所有变量及其属性的区域被各相共享并且代表了体积平均值，只要每一相的体积分数在每一位置是可知的，就能得出不同位置不同时刻下各相的分布情况。在结构物穿过水和空气的交界面时，涉及到空气和水两相流动的问题。在每个单元中，用a1代表空气的体积分数，a2代表水的体积分数，通过求解任何一相体积分数的连续方程即可实现对水和空气之间界面的追踪，知道不同时刻下液面的变形情况。

笔者基于VOF模型，在Fluent中运用动网格技术，将控制回转体运动的自定义函数(UDF)与流动控制方程耦合，模拟得到回转体在出水过程中的压力场和速度场。然后与实验及文献比较，在物体头部带水和尾部拖水，涡量分布方面都吻合得较好，为出水过程的进一步研究奠定了基础。

1 数值计算方法

1.1 计算模型及边界条件

为了便于与实验进行比较，笔者的计算模型类似张军的TR-PIV实验，只是运动速度改为5.1 m/s恒定不变，比其扩大了10倍，但这不影响压力、速度等场量的分布情况。

计算模型是一钝头回转体，头部直径D=30 mm，高 h=60 mm，计算域为 500 mm×800 mm，水深500 mm。网格划分在 Fluent专业前处理器Gambit上完成。网格采用非结构化三角形网格，在回转体周围进行网格加密，网格总数为81 822个，计算时间为17.5 h，如图1所示。

边界条件设定:回转体、计算域两侧和底部为固壁面无滑移条件，上端开口为压力出口条件，通过外加UDF控制回转体匀速直线运动。

1.2 控制方程

结构物出水是一种非定常绕流问题，控制方程主要有连续性方程和N-S方程。

连续性方程为:

N-S 方程为[10]:

由于回转体在出水过程中可能导致负的正压力，为了更好地模拟湍流的物理规律，选择Realizable κ-ε模型，关于κ和ε的输运方程如下[11]:

其中，σκ=1.0，σε=1.2，C2=1.9，由 Fluent湍流模型在控制面板上提供。C1=max[0.43，η/

VOF模型通过求解水或空气体积分数的连续方程来实现相与相间的界面追踪，对第q相，其方程如下:

式中，mqp为第q相对第p相的质量传输;mpq为第p相对第q相的质量传输;Saq为源项。主相容积比率的计算基于如下约束:

如果水和空气分别用下标1、2表示，则每一单元的密度由式(7)给出:

1.3 计算参数设置

控制动网格变化采用弹簧近似光顺模型和局部重划模型，并激活尺寸函数控制，防止计算过程中出现扭曲度过大的网格。压力和速度耦合采用PISO算法。壁面处理采用增强的壁面函数法。

2 计算结果和分析

图2为钝头回转体出水过程的速度矢量场，按等时间间隔选取了6个不同时刻。从图2中可以看出，在开始阶段，速度的分布范围很广，囊括了以回转体最大尺度为直径的圆，头部由于排挤作用，速度呈放射状，最大速度出现在底部中心区域，左右两侧速度较小，并有明显的回流区，左侧为逆时针，右侧为顺时针，这是回转体对水的拖带作用结果。随着物体的不断上升，尾部速度分布范围拉长，两侧的回流区也越来越不明显，在物体穿越水面进入空气后(t=0.082 s以后)，沿两侧壁面出现较大速度，是因为回转体将水带出，水受黏性及惯性力作用继续向上运动。

图3为钝头回转体出水过程的涡量分布图，同速度矢量场一样，选取了6个时刻，观察其变化过程。从图3中可以看出，在回转体底部两个拐角处分别分布一涡，并随着运动继续，涡变为带状，并逐渐从壁面脱落，形成涡列。在头部也可以看到左右分布着两个涡。

回转体在出水的瞬间，流体会因为黏性及惯性作用，附着在壁面上随物体一起运动，形成薄薄的一层水膜，之后由于侧面压力减小，水受到的粘压力不足以弥补重力，流体会慢慢滑落、凹陷，从而水膜破灭，壁面与空气接触。图4为回转体出水瞬间及以后的两相图。从图4中可以看出回转体拖带水的情况及变化趋势，开始阶段，头部带出大量液体，呈山体状，随着物体继续上升，水体逐渐变细，中部向内凹，最后呈带状。水膜破灭是从两侧壁面开始，先是出现一些不规整的小突起，然后逐渐扩大到头部，头部的水分又由于重力的作用滑落到侧壁，最后同底部带状水体融合。

图2 不同时刻下回转体的速度矢量图

图3 不同时刻下回转体的涡量图

图4 回转体出水瞬间及以后的两相变化示意图

回转体在整个出水过程中，会同时受到压力和流体黏性力的作用。计算过程中发现，两侧的黏性系数及底面的压力系数一直很小，几乎为零，而回转体在出水过程中所受总的流体力变化及顶部的压力系数如图5、图6所示。从图中可以看出，总流体力和顶部压力系数变化趋势一致，在水中起伏很大，当回转体穿过水面进入空气后，两者趋于稳定，在很小范围内浮动。

同张军等的TR-PIV实验结果比较可以看出，钝头回转体在出水过程中所引起的速度场、涡量场分布是吻合的，出水瞬间头部、尾部拖带水分情况也符合实际情况，这说明计算结果是准确的，可以为出水过程的进一步研究提供依据。

图5 回转体在出水过程中所受总流体力变化情况

图6 回转体在出水过程中顶部压力系数变化情况

3 结论与展望

笔者利用VOF模型，通过编写自定义函数(UDF)控制回转体运动，计算得到了钝头回转体在出水过程中的速度矢量场、涡量场的分布情况，以及出水后物体所带出水分的变化情况，为结构物出水的深入研究提供了依据。

结构物出水是一个涉及到两相流、流固耦合的非线性复杂问题，有很多情况有待深入研究。从涡量分布结果看，开始阶段及往后较长时间内，底部两个拐角处一直分布着一对点涡，这是否与物体形状、流体性质等有关还有待研究。

[1] 杨晓光.潜射导弹水下发射及出水过程三维数值研究[D].哈尔滨:哈尔滨工业大学图书馆，2009.

[2] KARMAN V T，WATTENDORF F L.The impact on seaplane floats during landing[J].National Advisory Committee for Aeronautics，1929(TN321):353-358.

[3] HAVELOCK T.The forces on a circular cylinder submerged in a uniform stream[J].Proc R Soc London Ser A，1936(157):526-534.

[4] LIJU P Y，MACHANE R，CARTELLIER A.Surge effect during the water exit of an axis-symmetric body traveling normal to a plane interface:experiments and BEM simulation[J].Experiments in Fluids，2001(31):241-248.

[5] 鲁传敬.轴对称细长体垂直出入水[J].水动力学研究与进展，1990，5(4):387-393.

[6] 蔡军辉，何友声，叶取源.平头圆柱体垂直出水粘性流动的数值模拟[J].水动力学研究与进展，1990，5(4):115-121.

[7] 张军，李英浩，金朋寿.垂直及斜出水流场的二维及三维 TR-PIV 试验[J].船舶力学，2005，9(2):9-17.

[8] 童亮，于罡，彭政，等.基于VOF模型与动网格技术的两相流耦合模拟[J].武汉理工大学学报:信息与管理工程版，2008，30(4):525-528.

[9] 江帆，黄鹏.Fluent高级应用与实例分析[M].北京:清华大学出版社，2008:140-164.

[10] 夏国泽，马乾初.船舶流体力学[M].武汉:华中科技大学出版社，2003:42-93.

[11] 王福军.计算流体动力学分析:CFD软件原理与应用[M].北京:清华大学出版社，2004:105-126.