﻿ 基于VOF法的不同形状液舱晃荡数值模拟
 舰船科学技术  2022, Vol. 44 Issue (8): 84-89    DOI: 10.3404/j.issn.1672-7649.2022.08.017 PDF

Numerical simulation of sloshing in tanks with different shapes based on VOF method
WU Yi-wei, ZHOU Jun, LI Hui, BI Zhou
Kunming Branch of the 705 Research Institute of CSSC, Kunming 650032, China
Abstract: Taking the liquid sloshing problem in the horizontal tank of an underwater vehicle as the research object, a two-dimensional horizontal tank numerical model is established based on the VOF method to study the response of different external excitations to the tank sloshing. The sloshing state of the free liquid surface under two special working conditions is calculated by simulation, and the errors of numerical simulation and literature results are analyzed, which proves the accuracy of the model calculation. Based on the calculation of the rectangular tank, the simulation analysis of the semi-circular-rectangular tank similar to the actual shape of ocean engineering is further carried out. According to calculation , Analyzing the influence of each influencing parameter on tank sloshing, the wave height duration curve in the ballast tank and the pressure time curve at a monitoring point on the bulkhead are obtained. There exist small numerical simulation’s error, it provides a theoretical basis for reducing the impact of sloshing on the underwater vehicle.
Key words: semicircle rectangular tank     liquid sloshing     STAR-CCM+     VOF method
0 引　言

Abramson[1]分析了多种不同类型液舱内的液体运动，研究各个参数对液面晃荡造成的响应。Faltinsen[2]构造了自由液面在矩形液舱内液体遭受外部（主要为余弦）激励的线性解析解，并拓展为同形状容器中自由液面晃荡问题的非线性解析解[3]，其研究成果被Faltinsen等[4]深层次应用于数值模拟更加典型的非线性计算。Cho等[5]模拟了受外部水平激励时矩形液舱的剧烈振荡。Wang等[6]对随机激励下液面的非线性水动力系数进行了计算。Liu[7]选取创新的数值模型去研究计算存在不完整自由液面的非线性晃荡。Lu[8]耗费了较长时间模拟液舱晃荡，从而得到较为收敛的数值计算结果，并简单阐述了能量耗散的原理。近年来数值模拟[9]的方法在不断升级进步，也出现了一些无网格化的方法，但采用大型商用软件针对晃荡问题进行的研究不多，为研究大型商用软件对晃荡问题的模拟能力及可靠性，本文采用STAR-CCM+软件建立2种不同形状的二维水平液舱模型，数值模拟外部余弦激励下的自由液面晃荡，其结果对评估水下航行器的稳定性具有借鉴意义。

1 晃荡问题的数值方法 1.1 条件假设

1）航行器稳定直航；

2）液面不可压缩，且为黏性流；

3）液舱不在冲击载荷作用下发生变形。

1.2 VOF法基本原理

 $\frac{{\partial \rho {u_i}}}{{\partial {x_i}}} = 0 ，$ (1)
 $\frac{{\partial \rho {u_i}}}{{\partial t}} + \frac{{\partial \rho \left( {{u_j} - u_j^m} \right){u_i}}}{{\partial {x_j}}} = \rho {f_i} - \frac{{\partial p}}{{\partial {x_i}}} + \rho v\frac{\partial }{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)。$ (2)

VOF(volume of fluid)方法是一种典型的界面追踪方法，它建立在欧拉网格的基础之上。该方法选取流体体积分数为界面函数 $F$ （在多相流中，称色函数，用 $C$ 表示） 这一变量来实现对计算域内相间界面的追踪。

$F$ 满足方程 $\dfrac{{{\rm{d}}F}}{{{\rm{d}}t}} = 0$ ，结合连续方程，则 $F$ 应满足

 $\frac{\partial }{{\partial t}}\left( {\theta F} \right) + \frac{\partial }{{\partial x}}\left( {\theta uF} \right) + \frac{\partial }{{\partial y}}\left( {\theta vF} \right) = 0。$ (3)
2 液舱晃荡数值模型建立和验证

2.1 仿真流程和计算参数设置

 图 1 仿真流程 Fig. 1 Simulation process

2.2 二维水平液舱模型及液体晃荡固有频率分析 2.2.1 二维水平液舱模型

 图 2 半圆-矩形STAR-CCM+网格示意 Fig. 2 Diagram of STAR-CCM+ grid of semicircle rectangular tank
2.2.2 液体晃荡固有频率分析

 ${\omega _i} = \sqrt {g{k_n}\tanh {k_n}d},\;\; n = 0,1,2 \cdot \cdot \cdot ,{k_n} = (2n + 1){\text{π}} /L 。$ (4)

Faltinsen[2]构造了势流解析解，给定水平液舱遭受一个外部激励： ${u_e} = - A\cos \omega t$

 $\begin{split} \eta =& \frac{1}{g}\sum\limits_{n = 0}^\infty {\sin \left[ {\frac{{(2n + 1){\text{π}} }}{{2a}}x} \right]} \cosh \left[ {\frac{{(2n + 1){\text{π}} }}{{2a}}h} \right]\times\\ &\left( { - {A_n}{\omega _n}\sin {\omega _n} - {C_n}\omega \sin \omega t} \right)- \frac{1}{g}A\omega x\sin \omega t，\end{split}$ (5)
 $\omega _n^2 = g\dfrac{{(2n + 1){\text{π}} }}{{2a}}\tanh \left[ {\dfrac{{(2n + 1){\text{π}} }}{{2a}}h} \right]，$ (6)
 ${A_n} = - {C_n} - \dfrac{{{K_n}}}{\omega }，{C_n} = \dfrac{{\omega {K_n}}}{{\omega _n^2 - {\omega ^2}}}，$ (7)
 ${K_n} = \frac{{\omega A}}{{\cosh \left[ {\dfrac{{(2n + 1){\text{π}} }}{{2a}}h} \right]}}\frac{2}{a}{\left[ {\dfrac{{2a}}{{(2n + 1){\text{π}}}}} \right]^2}{( - 1)^n}。$ (8)

2.2.3 数值模拟可靠性验证

 图 3 b=0.01 m， $\omega$ =0.5 $\; {\omega _0}$ 时舱壁处波高历时曲线 Fig. 3 Wave height time progress curve at bulkhead when b=0.01 m， ${\omega _0}$ ＝0.5 ${\omega _0}$

 图 4 b=0.000 4 m， $\omega$ =0.95 ${\omega _0}$ 时舱壁处波高历时曲线 Fig. 4 Wave height time progress curve at bulkhead when b=0.0004 m， ${\omega}$ ＝0.95 ${\omega _0}$

 图 5 自由液面初始云图 Fig. 5 Initial cloud diagram of free-surface
3 晃荡数值模拟 3.1 振幅工况对晃荡的影响

 图 6 矩形液舱模型及其网格划分 Fig. 6 Rectangular tank model and its mesh generation.

 图 7 矩形液舱不同振幅自由液面相对波高分布 Fig. 7 Relative wave height distribution of free-surface at different amplitudes under rectangular tank.

 图 8 半圆-矩形液舱不同振幅自由液面相对波高分布 Fig. 8 Relative wave height distribution of free-surface at different amplitudes under semicircle rectangular tank

$\omega$ =0.5 ${\omega _0}$ 工况振幅为0.003 m和0.004 m时， $\eta /b$ 值分别稳定在0.34和0.36； $\omega = {\omega _0}$ 工况振幅为0.003 m和0.004 m时， $\eta /b$ 值分别稳定在6.58和5.67。在 $\omega$ =0.5 ${\omega _0}$ 工况下，随着振幅的增加 $\eta /b$ 略微增加，增幅为5.9%；在 $\omega$ = ${\omega _0}$ 工况下，随着振幅的增加 $\eta /b$ 略微减小，减幅为13.8%。

 图 9 不同激励振幅下压力时间历程曲线 Fig. 9 Pressure-time history curve under different excitation amplitudes
3.2 频率工况对晃荡的影响

 图 10 不同频率自由液面相对波高历时曲线 Fig. 10 Relative wave height duration curve of free-surface at different frequencies

 图 11 不同频率自由液面相对波高分布 Fig. 11 Relative wave height distribution of free-surface at different frequencies

 图 12 不同激励频率下压力-时间历程曲线 Fig. 12 Pressure-time history curve under different excitation frequencies
3.3 网格无关性验证

 图 13 加疏与加密后的网格示意图 Fig. 13 Mesh diagram after encryption and thinning

 图 14 不同网格数量自由液面相对波高分布 Fig. 14 Relative wave height distribution of free-surface with different number of mesh
4 结　论

1）对于本文研究的二维水平液舱的晃荡问题，计算2种特殊工况下自由液面晃荡状态，并将数值模拟计算和解析解结果相比较，2条曲线的误差较小，证明了本文计算选取的液舱模型的可靠性。

2）矩形和半圆-矩形液舱模型下， $\omega = {\omega _0}$ 时，波峰的最小值随着激励振幅增大而增大，最大值减小，曲线总体呈现向上偏离的趋势； $\omega$ =0.5 ${\omega _0}$ 时，随着数值模拟计算的推进，曲线的最大相对波高值在下降一段时间后收敛。2种形状液舱模型均呈现了一段时间的非线性现象，随后幅值趋于稳定。不同振幅的曲线差别较小，这证明该情况下施加的激励晃荡是与线性相接近的。

3）通过对比矩形液舱模型和半圆-矩形液舱模型不同振幅工况下的自由液面波动规律，半圆-矩形液舱模型波高波峰波谷差值略大，波动周期略长。总体而言，2种液舱模型晃荡时自由液面波高总体变化规律较为接近，吻合度较高。

4）半圆-矩形液舱模型下，在位移振幅工况选取0.003 m和0.004 m的情况下，舱内液体晃荡的冲击会对舱壁施加压力，且此压力有周期性规律；在开始的40 s内，均呈现了较为显著的非线性现象，随后压力幅值趋于一个稳定值。且随着激励振幅的增大，舱壁监测点受到的压力波峰值减小而波谷值增大，波峰波谷的差值逐渐稳定。在0.003振幅工况下，随着激励频率的增大，液舱晃荡对舱壁上监测点压力的稳定范围减小。

5）半圆-矩形液舱模型下，在不同的激励频率工况下，在固有频率附近自由液面波高值达到最大值；激励频率小于或大于固有频率时， $\eta /b$ 值随着 $\omega$ 远离固有频率值而减小。在60 s后波高曲线趋于收敛，其波高历时曲线与激励曲线运动规律一致，呈简谐波动规律。在给定的频率工况达到水舱中流体的固有频率工况时，就会激发出最为强烈的液体晃荡，其晃荡幅值远远超过其他频率工况。因此在海洋工程研究过程中，应尽可能地控制好外部施加的频率使其远离固有频率，才能保障航行器航行时的安全。

5 结　语

 [1] ABRAMSON H N. The dynamic behavior of liquids in moving containers[J]. National Aeronautics and Space Administration, 1966. [2] FALTINSEN O M. A numerical nonlinear method of sloshing in tanks with two-dimensional flow[J]. Journal of Ship Research, 1978, 22: 193-202. DOI:10.5957/jsr.1978.22.3.193 [3] FALTINSEN O M, ROGNEBAKKE O F, LUKOVSKY I A, et al. Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with finite water depth[J]. Journal of Fluid Mechanics, 2000, 407: 201-234. DOI:10.1017/S0022112099007569 [4] FALTINSEN O M, TIMOKHA A N. Adaptive multimodal approach to nonlinear sloshing in a rectangular tank[J]. Journal of Fluid Mechanics, 2001, 432: 167-200. DOI:10.1017/S0022112000003311 [5] CHO J R, LEE H W. Numerical study on liquid sloshing in baffled tank by nonlinear finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(23/24/25/26): 2581-2598. [6] WANG C Z, KHOO B C. Finite element analysis of two-dimensional nonlinear sloshing problems in random excitations[J]. Ocean Engineering, 2005, 32: 107-133. DOI:10.1016/j.oceaneng.2004.08.001 [7] LIU D, LIN P Z. A numerical study of three-dimensional liquid sloshing in tanks[J]. Journal of Computational Physics, 2007, 227(8): 3921-3939. [8] LU L, JIANG S C, ZHAO M, et al. Two-dimensional viscous numerical simulation of liquid sloshing in rectangular tank with/without baffles and comparison with potential flow solutions[J]. Ocean Engineering, 2015, 108: 662-667. DOI:10.1016/j.oceaneng.2015.08.060 [9] WU G X, MA Q W. Numerical simulation of sloshing waves in a 3D tank based on a finite element method[J]. Applied Ocean Research, 1998(20): 337-355. [10] 郭海宇, 张志国, 冯大奎, 等. 基于VOF法的半球-圆柱体液舱晃荡数值模拟[J]. 中国水运, 2014, 14(2): 80-83.