基于LES和BPBE的舰船气泡尾流数值分析方法

金良安，闫雪飞，迟 卫，陈明荣

（海军大连舰艇学院 航海系，辽宁 大连 116018）

0 引 言

舰船尾流具有存活时间长，目标特征明显和易被探测等特点，因此一直备受关注，是军事领域的一个重要研究热点[1]。实验探测、理论计算和数值模拟是当前对尾流研究的三种主要方法，其中数值模拟由于具有无需设备、经济适用和节省人力物力等诸多优点而倍受青睐。

随着计算机性能的飞速提升，对尾流的数值模拟方法也在不断地进步和发展。舰船尾流中可供探测的特征参数包括速度、湍动能、耗散率及气泡尾流等，其中主要以气泡尾流的研究为主。田恒斗[2]通过建立舰船远程尾流场中单个气泡上浮运动与传质速率的数学耦合模型，分析研究了气泡尾流的存活时间、气泡数密度的衰减规律以及不同尺寸气泡相对数密度的分布特征。Stewart[3]采用大涡模拟技术（LES）和拉格朗日粒子跟踪方法对尾流环境下的气泡动力学进行了研究，重点分析了气体传质与扩散对气泡群尺寸演变的影响。朱东华[4]在舰船尾流相关特征参数求解的基础上对尾流的气泡数密度分布进行了三维数值模拟，该方法基于波尔兹曼型气泡输运方程和雷诺平均方程，计算出了尾流中不同尺寸的气泡数密度分布特征，发现尾流中的气泡尺寸主要以10～200 μm为主。其它针对舰船尾流气泡分布特征的数值分析研究更是屡见不鲜[5-9]。

当前针对舰船尾流的数值计算方法主要以直接数值模拟（DNS）和求解复杂的气液两相流模型为主，要消耗巨大的存储空间和计算时间，不是特定性能的巨型机根本无法实现，更不能达到实时性监测的军事斗争需求，为此，本文首次利用LES和BPBE方程对流场速度和气泡运动分开求解，大大地节省了时空复杂度，实现了舰船气泡尾流的大范围数值模拟（模拟长度为3 000 m）。

1 舰船尾流LES数值模型

1.1 模型建立的需求

针对当前气泡尾流主要数值研究方法的复杂性以及要占用大量的运行时间和存储空间的缺陷，将流体参量与气泡运动分开再来进行求解，不仅可以满足误差要求，还可节省大量的时空复杂度，具有事半功倍的效果。且流体计算和气泡运动方程都分别具有较好的移植性，便于用于其它两相流领域的应用和求解。为此，先建立尾流LES控制方程，以计算BPBE方程所需的一组速度参数，是气泡尾流求解的基础。

1.2 模型建立的过程

首先建立尾流模型的三维直角坐标系，沿尾流速度方向为x轴，纵向为y轴，垂直水面向下为z轴，对应各坐标轴的速度分量分别为u、v、w。使用格子滤波函数对N-S方程进行滤波，滤波后得到完整的LES控制方程如下式（其中i，j=x，y，z，展开方式满足求和约定）：

式中：ν是运动粘度系数，等式右边出现不封闭项τij，即为亚格子尺度应力，需要构建模型予以封闭，常用的Smagorinsky涡粘模型如下所示：

式中：cs为Smagorinsky常数，其值在0.1～0.27之间，Δ是格子滤波的长度量纲。在给出初始边界条件的情况下，采用数值离散方法即可对上述LES控制方程实现求解。

2 BPBE数值模型

2.1 模型构建的基础

BPBE方程通常用于鼓泡塔等湍流流场的气泡群运动特征的求解，因此可以移植过来应用于近尾流湍流流场的气泡群求解，而且包含的影响因素较为全面，下面将给出尾流中BPBE方程的详细构建过程。

一个典型的BPBE方程如下式所示[10]：

式中：n（z,d,t）为与空间 （z）、气泡尺寸 （d）、时间 （t）有关的单位体积气泡数密度函数，u是对应的 （z）方向的气泡速度，d（z,d）为与空间（z）、尺寸（d）有关的表示气泡尺寸因气体传质而引起变化的函数，此外，与气泡聚并有关的源项S（z,d,t）可以如下表示：

其中，

分别表示体积小于Vi的气泡与体积为（Vi- Vj）的气泡聚并为体积为Vi的气泡速率，气泡Vi与其它气泡聚并而导致气泡Vi消失的速率。这里假设聚并过程没有体积亏损[10]。

下面我们将基于舰船尾流实际情形对方程（3）进行适当的化简和离散处理，以方便程序更好更快地计算和求解。

2.2 模型建立的相关求解

仅考虑气泡数密度在轴向x和竖直方向y的变化，空间项可以如下表示：

注意到以下两个关系式：

则有：

对于方程（3）中的扩散与传质，有：

由（9）式可知，时间项可以被消去，因此（11）式可以化简为如下形式：

其中：ψ（d）表征扩散与传质，是关于直径d的一元函数。联立各式且考虑到气泡尺寸小到可以忽略分裂项，最终得到如下模型：

碰撞概率θij由下式给出[12]：

聚并概率由下式给出[12]：

其中： 初始膜厚度 h0＝5×10-4， 破裂液膜厚度 hf=1×10-8，d 为气泡直径，ρg为气体密度，ρl为液态密度，u是运动粘度，ξ为湍动能耗散率，Patm为标准大气压，CA为海水中的气体质量浓度，CI为气液界面处的气体质量浓度，具体公式见参考文献[1]。n（x,d）表示x处尺寸为d的单位体积气泡数密度。u（x,d）表示x方向尺寸为d的气泡速度。气泡群平衡方程最终化为只与坐标x和气泡尺寸d有关的二维偏微分方程，在求解前需要对轴坐标x和气泡尺寸d进行离散化处理。

3 计算与讨论

首先采用前文提出的LES模拟方法对流场速度进行求解，在已知初始条件的情况下即可实现对流场速度的数值模拟。流场初始流动参量的方法可采用MINER[13]提出的方法，湍动能耗散率的解算方法将结合MINER与LES计算结果进行求解。得到速度场后（展向、垂向、流向尺寸为：5 m×1 m×3 000 m），将其代入到BPBE方程求解BND的分布。本文将以文献[4]给出的美国某型驱逐舰为实船模型，各参数具体如下：舰长为155 m，舰宽为20 m，吃水6 m，排水量8 000 t，航速15 kns，螺旋桨直径5 m。

图（1）是x-y平面的位于流向位置（x=20 m）处的速度u等值线图，从图中可以看出速度在对称面处的值最大，沿展向和垂向逐渐减小，且呈对称分布。同时由图形的形状可以看出其分布与Kelvin尾流横截面结构相似。图（2）反映了某一位置处（y=0 m、z=1 m）速度u沿流向的变化，可以看出尾流平均速度非常小，在500 m以内以光滑的指数曲线迅速衰减，之后变化缓慢，与文献[14]测量结果的变化规律基本一致，这表明了本文LES计算方法的准确性，虽然在大小上有一定的偏差，造成误差的主要原因是算例中选取的目标船与文献[12]的实测目标船2者船型结构参数的不同引起的。

图1 速度u的等值面图Fig．1 The isoline-section of velocity u

图2 速度u沿流向x的变化Fig．2 The change of velocity u with x

3.2 尾流场BND计算与分析

为了进行数值计算，气泡尺寸需要按照大小进行分组，根据理论分析和实测经验，分组方法将如表1所示，分别给出了各组气泡的尺寸以及气泡数密度及单位体积分数，这样共得到15个相互耦合的差分方程组。模型将在普通台式机上运行，运行时间不超过5 s，相对以往文献的数值模拟方法大大提高了计算效率，可以达到实时监测的要求。

dbi /μm ni (1/m3) εdi (%) dbi /μm ni (1/m3) εdi (%) dbi /μm ni (1/m3) εdi (%)10 21 35 53 77 2.61e+10 6.31e+9 2.97e+9 1.69e+9 1.04e+9 0.010 9 0.024 5 0.052 4 0.106 0.204 109.6 151 203.5 268 346.5 6.01e+8 2.95e+8 1.51e+8 8.08e+7 4.5e+7 0.332 0.426 0.531 0.651 0.784 440.5 551.5 680.5 829.5 1 000 2.6e+7 1.56e+7 9.56e+6 6.06e+6 5.13e+6 0.931 1.09 1.26 1.45 2.14

本文重点研究气泡群BND沿流向的变化及相对分布的特征，取LES计算结果具有代表性的某一展向位置处（y=0 m）、某一垂向深度处（z=1 m）的一组流向速度参量作为气泡群平衡方程的流场速度。

图3 气泡数密度沿x方向的变化Fig.3 The change of bubble BND with x

图4 不同位置处气泡相对数密度Fig.4 Relative BND at different location

图（3）是BND总量沿x方向的分布，可以看出BND在离船尾500 m距离以内下降速度很快，之后衰减速度趋于缓慢，总体衰减规律与幂乘曲线近似。图中黑方块为文献[4]中的实测数据，可以看出两者形状基本相同。图（4）反映了不同位置处（x=500 m、x=1 500 m、x=3 000 m）的不同尺寸气泡相对BND的分布情况。由图可知尾流中的气泡成分主要以尺寸小于200 μm为主。在离船尾越近的地方，小气泡的相对BND越大，随着距离的增加，峰值开始向200 μm靠拢，但这种变化不大。在x=3 000 m的位置处，BND最大的气泡尺寸约为70 μm。这些结论都与相关文献和实验测量的结论[4-5]比较接近，进一步表明了本文使用BPBE方程对舰船气泡尾流进行求解的可行性。

4 结 语

舰船尾流一直是国内外研究的热点，尤其是气泡尾流的研究，在海洋工程、舰船安全、雷达监测以及武器制导等诸多领域，都具有重要的理论应用价值。针对当前对舰船尾流进行数值模拟的现有方法的复杂性及计算效率低下的缺陷，本文基于LES模拟和修正的BPBE方程，实现了尾流速度和气泡运动的分开求解，大大提高了计算效率，且具有普遍适用性，在普通台式机上即可实现模型的求解。计算结果表明：BND沿尾流流向的分布基本符合指数分布，在500 m以内的区域衰减速度很快，500 m以外的区域衰减速度缓慢；在离船尾越近的地方，小气泡的相对BND越大；在3 000 m距离处BND最大的气泡尺寸约为70 μm。所得结论与相关文献比较接近，可望为舰船气泡尾流特征研究提供新的有效手段。

[1]田恒斗,金良安,迟 卫.船舶远程尾流场中气泡上浮运动模型[J].船舶力学,2010,14(11):1195-1201.Tian Hengdou,Jin Liangan,Chi Wei.Rising process model of bubble in the far wake of ship[J].Journal of Ship Mechanics,2010,14(11):1195-1201.

[2]田恒斗,金良安,王 涌,等.考虑单气泡运动特性的舰船尾流气泡分布研究[J].兵工学报,2011,32(9):1126-1130.Tian Hengdou,Jin Liangan,Wang Yong,et al.Study on wake bubble distribution based on single bubble motion characteristics[J].Acta Armamentarii,2011,32(9):1126-1130.

[3]Stewart M B,Miner E W.Bubble dynamics in a turbulent ship wake[R].NRL Memorandmm Report,6055,1987.

[4]朱东华,张晓晖,顾建农,等.舰船尾流及其气泡数密度分布的数值计算[J].兵工学报,2011,32(3):315-320.Zhu Donghua,Zhang Xiaohui,Gu Jiannong,et al.Numerical calculation of ship wake and its bubble number density distribution[J].Acta Armamentarii,2011,32(3):315-320.

[5]顾建农,田雪冰,张志宏,等.舰船主尺度参数对舰船气泡尾流几何特性的影响[J].舰船科学技术,2012,34(1):6-10.Gu Jingnong,Tian Bingxue,Zhang Zhihong,et al.The effect of the ship main scale parameter upon ship bubble wake geometry characters[J].Ship Science and Technology,2012,34(1):6-10.

[6]陈圣涛,王慧丽,王运鹰,等.舰船气泡尾流特性的数值模拟和实验研究[J].船舶力学,2012,16(4):342-349.Chen Shengtao,Wang Huili,Wang Yunying,et al.Numerical simulation and experimental study of a ship’s bubble wake[J].Journal of Ship Mechanics,2012,16(4):342-349.

[7]杨 力,刘慧开,等.舰船尾流气泡特征的运动规律研究[R].武汉:海军工程大学,2005.Yang Li,Liu Huikai,et al.Movement regularities study of Bubbles in Ship Wakes[R].Wuhan:Naval University of Engineering,2005.

[8]Smirnov A,Celik I,Shi S.LES of bubble dynamics in wake flows[J].Computers&Fluids,2005(34):351-373.

[9]Carrica P M,Drew D,Bonetto F,et al.A polydisperse model for bubbly two-phase flow around a surface ship[J].International Journal of Multiphase Flow,1999,25:257-305.

[10]Ryszard P W,Moniuk W l,Bielski P,et al.Modelling of the coalescence/redispersion processes in bubble columns[J].Chemical Engineering Science,2001,56:6157-6164.

[11]徐麦容,刘成云.水中浮升气泡的半径和速度变化[J].大学物理,2008,27(11):14-17.Xu Mairong,Liu Chengyun.The change of radius and velocity of the rising bubble in water[J].College Physics,2008,27(11):14-17.

[12]Prince M J,Blanch H W.Bubble coalescence and break-up in air-sparged bubble columns[J].A.I.Ch.E.Journal,1990,36(10):1485-1499.

[13]Miner E W,Ramberg S E.Method for approximating the initial data plane for surface wake simulation[R].US:RL-MR-6376,1988.

[14]James K E,John R D,Brian A M.The radar image of the turbulent wake generated by a moving ship[J].London Research and Development,1993,755:343-346.

[15]石晟玮,蒋兴舟,石 敏,等.舰船尾流气泡运动特性研究[J].武汉理工大学学报（交通科学与工程版）,2007,31(5):764-767.Shi Shengwei,Jiang Xingzhou,Shi Min,et al.Movement characteristic of bubbles in ship wakes[J].Journal of Wuhan U-niversity of Technology(Transportation Science&Engineering),2007,31(5):764-767.