复杂流场下气泡界面不稳定现象的数值模拟方法综述

倪宝玉，张阿漫，张忠宇，李 帅

（哈尔滨工程大学 船舶工程学院，哈尔滨 150001）

复杂流场下气泡界面不稳定现象的数值模拟方法综述

倪宝玉，张阿漫，张忠宇，李 帅

（哈尔滨工程大学 船舶工程学院，哈尔滨 150001）

气泡界面不稳定现象一直是国内外关注和研究的热点。近年来，边界元方法（Boundary Element Method，BEM）被广泛应用于气泡动力学模拟，文章则侧重综述边界元法及其相关技术在模拟气泡界面不稳定现象中的应用与发展。首先，简要回顾不同边界不同背景流场下气泡溃灭坍塌的边界元法研究进展；其次，分别回顾和讨论气泡撕裂、融合和在自由面破裂三种典型界面不稳定现象的研究进展；再次，回顾和讨论了气泡模拟过程中可压缩性、表面张力和粘性三种因素的影响；最后，在上述基础上提出了一些尚需进一步解决的问题。该文旨在向国内同行学者介绍目前气泡界面不稳定现象的最新研究进展，希望能为国内学者研究气泡动力学提供一定参考。

气泡；界面不稳定；边界元法；撕裂；融合

1 引 言

气泡在流体中的运动学形态和动力学特性一直是流体力学、环境工程、医学和船舶与海洋工程等多种领域的研究热点之一，在自然界中随处可见流体中的气泡现象，在各工程领域的应用更不胜枚举。仅以船舶与海洋工程领域为例，较著名的气泡现象就有高速运转的螺旋桨周围形成的多种形态螺旋桨空泡[1]，有水下攻击武器爆轰后形成的水下爆炸气泡[2]，也有近年来绿色造船引入的气泡减阻技术中的微气泡或气穴[3]等。

气泡作为气液两相的自由交界面，在数学模型上属于非线性初值—边值问题。在数值模拟上，对于处理气泡这种带有自由交界面的问题，可以采用有限元法（FEM）或者边界元法（BEM）等方法。有限元法求解通常可分为两大类方法：界面追踪法（Interface Tracking Methods）[4]和界面捕捉法（Interface Capturing Methods）[5]，二者各有利弊。界面追踪法可以精确地定义流体界面，但是需要非常精密而复杂的网格重建工作，尤其当界面发生破损或融合等强非线性运动时。此外，计算过程中通常存在质量和能量损失。相反，以流体体积法（VOF）[6]、水平集法（Level-Set）[7]以及任意拉格朗日-欧拉法（ALE）[8]为代表的界面捕捉法在流域中弱化界面的概念，而是通过附加流体变量如质量分数、空隙度或比热率等来区分不同的流体。然而若想精确捕捉界面可能需要很多网格单元，同时在仿真过程中易出现界面模糊效应。

不同于有限元法，边界元法仅对流域的边界（如气泡表面）进行划分，而无需对整个流域内部及气泡内部进行划分，从而将问题的维度降低一维。BEM的最大优点在于高效性，网格和节点数的减少使BEM所消耗的时间相对于有限元法至少缩短几个量级，通常有限元法需要十几个小时或几天才能运行完成的结果应用边界元法几个小时就可完成。虽然边界元也有自己的不足，例如通常假设流体为势流，对于求解粘性有很大局限性，或无法计及物面内外的质量和能量交换等，但是由于其避免了重构交界面两侧流体网格的优势，在研究气泡这种瞬态自由面大变形问题方面，还是获得了巨大的应用和发展。

通常而言，气泡并不会稳定地存在于流体中，而是在流体环境或周围边界的复杂影响下表现出极不稳定现象，如剧烈溃灭、撕裂以及融合等。气泡界面的这种不稳定现象一方面会改变气泡自身的动力学行为，另一方面也会改变其对周围结构的载荷特性，故关于气泡界面不稳定现象的研究一直是国内外的热点之一。本文将回顾并讨论气泡不稳定现象的研究进展和技术难点，并主要侧重于边界元法（BEM）研究气泡不稳定现象的相关数值处理和技术手段。

2 气泡溃灭坍塌

基于势流理论的边界元法（BEM）或称边界积分法（BIM）在瞬态自由交界面上最早应用，可能是Longuet-Higgins和Cokelet（1976）[9]。他们应用第二格林函数将拉普拉斯方程转化为边界积分方程，并采用椭圆积分数值求解积分方程，从而成功地解决了非线性瞬态水波问题。此后，边界积分法迅速发展，在气泡动力学上应用广泛（Lenoir（1979）[10]，Guerri等人（1981）[11]，Blake 等人（1986，1987）[12-13]，Chahine（1990）[14]，Best（1993）[15]，Zhang 和 Duncan（1994）[16]，Brujan 等 人 （2002）[17]，Lee，Klaseboer 和Khoo（2007）[18]）。边界元法求解气泡动力学问题的基本思路如下：已知气泡表面的速度势，应用格林第二函数求解可获得气泡表面的法向速度∂φ/∂n，同时应用速度势差分法可获得表面切向速度∂φ/∂τ，进而可合成节点的全速度▽φ。在时间步进中，采用拉格朗日法追踪气泡表面节点，并应用拉格朗日坐标求解气泡表面伯努利方程，获得每一时间步内速度势增量Δφ。在获得新的速度势和新的气泡节点坐标后，即可进入循环，获得气泡随时间的推进演化历程。

在不同边界或不同背景流场条件下，气泡在收缩阶段可能会发生不稳定坍塌溃灭，并形成高速射流。最典型的两类边界是刚性壁面和自由液面。对于刚性壁面附近气泡的求解中，对于气泡表面，速度势已知，满足狄利克雷（Dirichelet）边界条件；对于刚性壁面，法向速度已知，满足纽曼（Neumann）条件，仅需对边界积分方程得到的系数矩阵重组即可获得所需未知量。研究发现气泡在壁面附近膨胀阶段会轻微远离壁面，而收缩阶段会被壁面强烈吸引，形成指向壁面的射流。壁面对气泡的这种效应称为“Bjerknes效应”（Bjerknes，1966[19]）。对于自由液面附近的气泡，自由面和气泡面均满足狄利克雷边界条件，在时间步进中也采用拉格朗日法，这样即可采用与求解气泡类似的方法追踪自由面的变形。自由面附近的气泡动态研究基本上是与壁面附近气泡同步展开的，然而与壁面不同，在不考虑浮力影响下，气泡在坍塌阶段将远离自由面，并形成指向背离自由面的射流，而自由面将形成明显的水冢现象（Blake和Gibson，1987[20]）。无论是刚性壁面还是自由液面，气泡在溃灭后期，射流都将穿透气泡另一侧壁面形成环状气泡。环状气泡属于双连通域，现在应用的较广泛的求解方法有“涡面模型”[17]和“涡环模型”[21]。关于边界附近气泡不稳定溃灭，尤其是环状气泡的研究方法和进展，请参见文献[22]和[23]。

对于不同背景流场，本文选择船舶与海洋工程中经常遭遇的两类背景流场为例，一是高速螺旋桨诱导的漩涡流场，一是水下爆炸引发的冲击波流场。

对于漩涡流场，经研究表明，螺旋桨诱导流场中气泡将随着涡流做螺旋线运动，在运动过程中自身还有可能撕裂为多个子气泡（Choi，2006[24]）。最早采用边界积分法研究漩涡流场中气泡的变形和运动的当属 Chahine（1983）[25]。Duraiswami和 Chahine（1992）[26]采用渐进展开法求解势流中某线涡附近气泡的变形，其中假设气泡本身的尺度相对于气泡与线涡的距离是个小量。数值计算表明气泡将比正常流中以更快的速度向下游运动，同时将形成指向涡心的斜射流。Hsiao和Pauley（1999）[27]应用雷诺平均NS方程模拟翼梢涡流场，并建立球状气泡模型研究空化初生现象。之后，Hsiao和Chahine（2003）[28]又进一步扩展球状气泡模型，研究Rankine线涡内空化气泡的尺度效应。然而，这种算法的缺点在于无法考虑气泡与流场的相互作用，也无法计及非球状效应。为了克服这个缺点，Hsiao和Chahine等人（2004，2006）[29-30]基于势流理论开发了非球状模型来研究空化初生、气泡坍塌溃灭和辐射噪声等，尤其是研究了气体扩散和溶解浓度的影响，获得了许多有益的结论。图1给出了文献[31]采用直接计算法（DNS）和文献[32]采用粒子追踪法（PTM）计算得到的螺旋桨流场中气泡的捕获过程。在国内，王献孚（2009）[33]基于修正的RP方程总结了螺旋桨流场中球状空泡初生判别方法。韩宝玉等人（2011）[34]应用Hsiao等人开发的球状气泡平均压力模型（SAP）预报了不同初始半径的气泡的空泡初生数，并预报了相应的辐射噪声。然而，国内直接应用边界积分法求解涡流场中气泡变形和撕裂等现象的文献尚十分缺乏。

图1 直接计算法（a）[31]和粒子追踪法（b）[32]计算气泡在漩涡流场的运动轨迹Fig.1 Bubble’s trajectory in a vortex field calculated by DNS(a)[31]and PTM(b)[32]

对于冲击波流场，除了重力场中流体的静压力梯度和边界效应（壁面或自由面）外，当气泡处于冲击波流场时，强大的冲击波压力也会促使气泡形成射流，且射流速度可能达到数千米每秒的量级（Bourne和Field，1992[35]）。Blake等人（1999）[36]应用边界元法建立了低频振动的刚性边界附近气泡的动态响应，他们在伯努利方程的重力项中附加了一项考虑了边界正弦振动引起的附加压力，此项附加压力是时间的函数但不是空间的函数，以此模拟不同振动频率和振幅强度下气泡的变形和射流等。Klaseboer等人（2006）[37]借鉴Blake等人的思想，在伯努利方程中引入随时空均变化的冲击参考压力，进而应用边界元研究了某一气泡与矩形脉冲的相互作用。在模拟中假设可压缩性仅在气泡坍塌后期影响明显，在大部分气泡运动周期内惯性项占主导作用。Klaseboer等人（2007）[38]又进一步考虑了某一碎石冲击波与气泡的相互作用，在数值模拟中考虑了球状气泡初始形态（如膨胀、静止或收缩）对于射流速度等影响。国内，倪宝玉等人（2011）[39]在Klaseboer等人的基础上，考虑爆炸冲击波在船底和水底形成的反射波对水下爆炸气泡的影响，将边界元理论计算应用到工程领域。图2即给出了文献[39]采用边界元法（BEM）与文献[40]采用自由拉格朗日法（FLM）数值模拟冲击波流场中气泡坍塌射流的演化过程。经过计算发现，反射波对气泡射流载荷的影响较弱，在工程计算中可将冲击波载荷和气泡射流载荷分别计算，无需考虑二者间的耦合作用。当水下武器多发攻击目标时，异相爆炸间也会存在冲击波与气泡载荷的相互作用，尚需进一步研究。

图2 自由拉格朗日法（上）[40]和边界元法（下）[39]计算冲击波下气泡坍塌对比图Fig.2 Contrast of bubble collapse under shockwave calculated by FLM(top)[40]and BEM(bottom)[39]

3 气泡撕裂和融合

除了上文阐述的溃灭现象和环状气泡现象，在日常生活中气泡不稳定界面中更为常见的有气泡的撕裂、融合，以及自由面处破裂等。本章则在上述的基础上进一步分析气泡撕裂和融合的动力学行为，尤其关注边界元法在处理撕裂和融合问题的处理手段。

3.1 撕裂

气泡撕裂是自然界中常见现象之一，导致气泡撕裂的真实物理原因可能十分复杂，这里总结以下几个方面：首先，较容易理解的稳态气泡遭受外力或外物干扰，Kim和Kim（1990）[41]应用不同截面形状的物体主动切割某一稳态上浮气泡，总结浮体的形状、速度等因素对气泡撕裂的影响；其次，较常见的是气体体积增大，当大于某一临界值后，气泡在微小扰动下即自发破裂，Wilkinson等人（1993）[42]研究了不同气体密度和不同液体特性对于气泡碎裂的影响；再次，在湍流漩涡作用下，由于粘性中心半径、速度环量和压力等在涡流场内变化剧烈，气泡十分容易撕裂，Martinez-Bazan等人（1999）[43]即研究了湍流场中气泡的运动与撕裂现象；最后，不同边界的作用也会导致气泡撕裂，例如在两个壁面的Bjerknes力的相互吸引下，一个气泡在收缩阶段即可能分裂为两个子气泡。

为了揭示气泡撕裂现象及其载荷特性，人们进行了大量的研究。在实验方面，Ishida等（2001）[44]应用水箱模拟气泡在两水平壁面间的对称撕裂与非对称撕裂，并通过压力传感器测得壁面与流场的压力变化；Choi和Chahine（2003）[45]应用高速摄影拍摄了气泡在两倾斜壁面和两垂直壁面间的撕裂过程。图3给出了采用电火花气泡模拟两变化距离刚性壁面内气泡的撕裂现象。

图3 两平行壁面间气泡撕裂现象实验Fig.3 Experiment on breakup of a bubble between two horizontal walls

数值模拟上，早期 Ellis（1956）[46]，Lauterborn 和 Hentschel（1985）[47]等人都曾研究过两刚性边界间单个气泡的运动，预测气泡收缩后期可能形成沙漏状；Mohammad（2006）[48]等系统地研究了轴对称气泡在两平行壁面间的颈缩状态，对气泡撕裂前的行为做了详细的工作，但没有模拟出气泡撕裂及以后的过程；Choi和Chahine等（2003）[45]提出了轴对称气泡的撕裂准则，建立了轴对称撕裂模型并模拟了对称撕裂与非对称撕裂；张阿漫和倪宝玉等人（2010）[49]在Choi和Chahine研究的基础上提出了三维气泡撕裂准则，研究了两刚性壁面间单个气泡的对称和非对称撕裂现象，给出了不同工况下气泡形态、射流速度、开尔文冲量等物理量的变化。

图4 气泡撕裂数值过程示意图（来自于文献[49]）Fig.4 Process of the split of a 3D bubble(data from Ref.[49])

不直接模拟气泡撕裂的真实过程，在数值上人工干预，重新分布满足分裂准则的气泡节点，使得一个气泡分裂为两个子气泡，图4给出了边界元法处理气泡撕裂的过程，具体实行过程如下：

（a）确定局部“对称轴”：如图4（a）所示，气泡在撕裂前期会发生较明显的“颈缩”现象，在局部呈现筒状结构，选择筒状结构的中心线为局部“对称轴”。

（b）计算最小距离；在每个时间步内，计算颈缩区节点到局部“对称轴”的距离，并将最小距离设定为Δs。当无量纲距离Δs≤0.03后，认为满足撕裂准则。否则继续下一步计算，直至满足撕裂准则为止，这个过程通过程序自动完成，无需人工干涉。

（c）确定“破裂线”；当Δs满足撕裂准则后，设定对应此Δs的节点编号为I1，与I1在同一圈的节点I2，…，Im则被定义为“破裂线”。“破裂线”上所有节点将全部移除，并将I1，I2，…，Im节点坐标的算术平均值作为新的一点J的空间坐标。

（d）撕裂中心点：将中心点J沿局部“对称轴”轴向上下各移动微小距离0.001，得到两个新节点J1和J2，并将原来断开的节点分别与J1和J2连接上，重新编号节点号和单元号，从而实现了一个气泡到两个子气泡的分裂，如图4（b）所示。

（e）计算新节点物理量；J1和J2点的速度势均取“破裂线”上点I1，I2，…，Im的算术平均值即可更新子气泡的速度势和表面位置。

3.2 融合

与撕裂相对，融合现象是气泡动力学研究的另一热点，在许多领域如气-液/泥浆反应堆、氧化/氢化反应、发酵、污水处理以及最近的细胞工程等均有广泛应用。通常而言，对于气泡融合可以化为三个阶段（Rungsiyaphornrat等人，2003[50]）：第一阶段，子气泡相互吸引，相对运动至较接近的位置，足以保证一段时间的演化后能达到融合条件；第二阶段，两个气泡在继续吸引下表面开始变平坦，其间的液体逐渐被挤压出去，很长一段时间内液膜越来越薄；当然，如果此时气泡在过大的碰撞力下分开而反向运动的话，则不会发生融合。第三阶段，当液膜变得很薄，进入某一临界值后，分子力等内力加之微小外界扰动则将促使薄膜破裂，气泡融为一体。Chen等人（2005）[51]给出了气泡融合三个阶段的示意图，如图5所示。

图5 气泡融合示意图（图形改自于文献[51]）Fig.5 The sketch of bubble coalescence(figure revised from reference[51])

与液滴一样，在两个气泡融合前的小段时间，其间形成的液膜厚度非常小，可达到微米数量级。处理液滴融合的薄膜问题时，人们常采用流体润滑理论模型（Chevaillier等人，2006[52]），这样可以很好地解决液体薄膜的计算。在模拟气泡融合时，起初也试图引入流体润滑理论模型，但经过对两个气泡的数值模拟，发现气泡融合前的液膜运动规律并不能用流体润滑理论控制。这是因为液膜的雷诺数[50]至少在100以上，有的甚至达到1万以上，远远不能满足流体润滑理论应用时雷诺数小于1的要求[53]。关于气泡融合的实验数据相对较少，Klaseboer等人（2000）[53]实验研究了匀速运动接近的两气泡的融合现象，而后，Rungsiyaphornrat等人（2003）[50]应用轴对称势流理论模型求解了两个水下爆炸气泡同相和异相融合，数值结果与Klaseboer等人的实验进行对比，吻合度良好，发现同相气泡融合后将形成相互吸引的对射流。张阿漫等人（2010）[54]在Rungsiyaphornrat等人的基础上，采用边界元法数值建立了同相气泡三维融合模型，数值模拟了两同相爆炸气泡水中融合现象，并对比融合气泡和单个气泡引起的滞后流和压力的差异。然而，上述模型均计算到对射流相互接触时刻为止，没有考虑对射流碰撞后的环状气泡和环状射流阶段。

图6 气泡融合数值过程示意图Fig.6 Process of the mergence of two 3D bubbles

与撕裂类似，并不真实模拟气泡融合过程，仅在数值计算过程中，当气泡网格节点距离满足融合准则后，进行人工干预，重新组合网格而形成一个新的气泡。图6给出了边界元法处理融合过程示意图，具体实行过程如下：

（a）计算最小距离：设气泡表面节点个数分别为N1和N2，则每个时间步内均计算两个气泡各节点间的距离，得到N1×N2个距离，并假定其中最小距离为Δs，对应的两个节点分别标注为I和I′，如图6（a）所示。当无量纲距离Δs≤0.03后，则认为满足融合准则，否则继续下一步计算，直到满足融合准则为止，这个过程通过程序自动完成，无需人工干涉。

（b）确定“融合线”：当获得了满足融合准则的节点I和I′后，以节点I为例，将与I相连接的第一圈节点I1，I2，…，Im定义为“融合线”。在实际计算中，如果网格尺寸较小，也可以将融合线扩展到I周围的第二圈第三圈节点。类似地，在I′附近也会找到对应的一条“融合线”。

（c） “融合线”缝合：将融合线上的节点 I1，I2，…，Im和 I1′，I2′，…，In′人为去除，并取对应节点坐标的算数平均值而得到新的一圈节点坐标J1，J2，…，Jk，其中k=max m,（）n 。假若m≠n，需在认为去除融合线上的节点之前，先自身插值，确保在取平均值时，所有节点均可寻找对应的节点。

（d）重新编号：分别将J1，J2，…，Jk与原来分离的节点连接上，并将节点号和单元号重新编号，至此两个子气泡就融合为一个气泡了，如图6（b）所示。

（e）计算新节点物理量；对于J1，J2，…，Jk点的速度势均取移除点的算术平均值即可更新子气泡的速度势和表面位置。

3.3 自由面处破裂

上浮气泡或接近自由面附近水下爆炸气泡均会产生在自由面的破裂现象。对于上浮气泡而言，实验和数值均发现气泡破裂后，会形成向上的小水滴，小水滴对气液两相间的质量交换起到重要作用，在气泡体积谱分析、通风、沸腾、去气、蒸馏以及蒸发方面均有重要应用。与气泡融合现象类似，通常认为气泡破裂的全部过程也可划分为三个阶段（Georgescu等人，2002[55]）：第一阶段，接近阶段。此阶段气泡加速上浮，越来越接近自由液面，同时自由面微微隆起。有趣的是，实验观察（Malysa等人，2005[56]）到气泡不是立即与自由面融合，而是在自由面附近微幅弹跳，直至振幅衰减，气泡在自由面下处于平衡位置；第二阶段，融合阶段。此阶段气泡上表面与自由面间形成的水膜逐渐变薄，直至厚度小于某个临界值，则水膜在微小扰动下将破碎；第三阶段：坍塌阶段。此阶段气泡内部将和大气连通，气泡将坍塌。此时，表面波会汇聚于气泡坍塌轴心并产生一股向上的高速射流。并形成指向上方的不稳定射流。图7即给出了气泡在自由面处破裂的示意图。在这三个阶段中，将可能形成两种类型的水滴[53]：其一发生在第二阶段水膜破碎时刻，破裂的水膜碎片将四处飞溅，形成所谓的“水膜液滴”；其二发生在坍塌阶段末期，此时指向上方的高速射流是不稳定的，在射流中上部可能发生连续的破碎，并产生若干“射流水滴”。射流水滴可能携带较大的能量，故学者们相对更关心射流水滴的变化。

图7 气泡在自由面处破裂示意图Fig.7 The sketch of bubble bursting at free surface

由于射流水滴可能携带较高的能量，本文主要关注“射流水滴”的相关进展。在实验上，Kientzler等人（1954）[57]即开始应用高速摄影仪观察气泡在自由面的破裂现象，并指出第一滴射流水滴的体积约为原始气泡体积的十分之一。Garner等人（1954）[58]描述了蒸馏水和酒精内射流水滴的速度和数量，指出原始气泡体积由几十微米增加到几毫米时，第一滴射流水滴的速度将由几十米每秒迅速降到几米每秒，同时形成的射流水滴的个数将减少，最多的射流水滴数目也不会超过10个。Resch和Afeti（1991，1992）[59-60]实验得到，对于海水而言，当原始气泡半径大于3 mm后，则无法形成射流水滴现象。在数值模拟上，由于自由面边界条件的强非线性，气泡在自由面处的破裂现象一直是数值建模方面的一大难点。最早的数值模拟研究应属Boulton-Stone和Blake（1993，1995）[61-62]，他们应用边界元法数值模拟气泡破裂过程，并引入边界层考虑弱粘性效应。研究表明，表面张力对于气泡在自由面处破裂起到十分重要的作用，而粘性对于前期气泡运动影响很弱，仅在射流形成后起到较小的耗散效应，阻碍射流水滴的形成。Sussman和Smereka（1997）[63]采用Level set法求解NS方程，数值模拟初始时刻轴对称球状气泡在水平自由面处破裂，但气泡半径均选取大于4 mm的工况，与之前的0.5-3 mm小气泡没有对比性。Georgescu等人（2002）[55]采用带有粘性效应修正的边界元法模拟气泡破裂过程的射流水滴现象，在时间步进上采用二阶泰勒展开。数值结果与实验值进行对比，表明此方法在模拟破裂前期精度很好，在射流水滴形成的后期误差有所增大，故得到结论是仅在整个气泡演化过程的后期NS方程才可能需要采用。

4 几类主要影响因素

气泡动力学进展中，学者们十分关心可压缩性、表面张力和粘性等对于气泡不稳定溃灭坍塌，撕裂融合以及相关载荷的影响规律，并试图揭示在什么样的条件下这些因素起主要作用必须考虑，或起微弱作用可以忽略。故在本节中，分别阐述国内外对这几类影响因素的研究进展。

4.1 可压缩性

Prosperetti和Lezzi（1986）[64]应用渐近匹配展开法对气泡壁马赫数展开到一阶项，从而研究在可压缩流体中气泡的径向运动。他们获得了描述弱可压缩流体内球状气泡运动的单一参数方程组，此方程组包含了著名的 Herring-Keller方程（Herring，1941[65]；Keller和 Kolodner，1956[66]）。 Prosperetti和 Lezzi（1987）[67]又进一步将马赫数展开到二阶项，并指出对于不可压缩流体的一阶项修正已经在很大程度上考虑了可压缩效应，二阶项仅有较弱的影响。Pozrikidis（2001）[68]采用边界积分法数值模拟某无粘二维可压缩气泡在Stokes流中的运动，发现在纯剪切流或纯拉伸流中，初始圆环状的气泡在收缩后最终均将达到一稳定态，且稳定态的面积比初始圆环状面积大。Darren和Crowdy（2003）[69]扩展了Pozrikidis的边界积分法，将原本的积分方程转化为两个耦合的非线性常微分方程组，此耦合方程组的计算更直接更稳定，且需要更小的计算机内存。数值模拟还发现绝热气泡较等温气泡更容易膨胀。Wang和Blake（2010，2011）[70-71]也采用渐近匹配法对气泡壁马赫数分别进行了一阶和二阶展开，研究可压缩流场中非球状气泡的运动。与Prosperetti和Lezzi不同的是，Wang和Blake将流场划分为内外两部分，对于接近气泡的内流域认为一阶和二阶均满足不可压缩的拉普拉斯方程，对于远离气泡的外流域直接采用一阶和二阶线性波动方程求解，在内外域交接面进行展开匹配。通过混合拉格朗日—欧拉法和修正边界元法求解气泡的变形，此方法的主要优点在于保证精度的同时提高了全局求解波动方程的效率。

Geers和Hunter（2002）[72]采用双渐近法研究了可压缩流场中水下爆炸气泡的球状运动，给出了气泡运动的经验公式。王诗平（2011）[73]在Geers和Hunter模型的基础上，采用双渐近法求解了可压缩流场非球状气泡的运动。Geers和Hunter采用频域内局部和全局近似，高低频匹配后再通过拉普拉斯变换到时域问题，而王诗平（2011）[73]采用时域内前期逼近和后期逼近，在时域内匹配求解，在数值实现上更为方便。在此基础上，王诗平研究了自由场和边界附近可压缩气泡的运动，并指出对于水下爆炸气泡第一周期的运动，可压缩性影响较小，工程应用中可以忽略。

4.2 表面张力

Beeching（1942）[74]首次在Rayleigh球状气泡运动中考虑了表面张力2σ/R的影响。此后，Akhatov等人（1997）[75]数值研究了表面张力和气体扩散对于声致发光气泡稳定性的影响，指出表面张力对于小气泡运动十分重要，可能是声致发光气泡稳定脉动的主要原因之一。在数值模拟上，表面张力对于气泡运动影响的难点在于局部曲率的求解，因为对于给定的流体和气体，表面张力系数是已知的。对于二维或者轴对称气泡而言，曲率有解析公式，可参见Bronshtein和Semendyayev（1997）[76]；对于三维气泡，曲率的数值求解是有一定难度的。Brackbill等人（1991）[77]开发连续谱法求解曲率，核心思想为将气液交界面认为是有限厚度的过渡区域，在过渡区域内定义与局部曲率成比例的力密度函数。连续谱法消除了对交界面重构的需要，简化了表面张力的计算方法，可精确求解二维或三维曲率，而无需受限于模型的网格数目、复杂度以及动态演化等因素。Zinchenko等人（1997）[78]阐述了一种更简便的求三角形单元表面曲率和法向量的迭代算法，该算法是在Rallison（1981）[79]思想上的进展，主要是通过最小二乘法拟合局部最优抛物面，数值计算表明该算法鲁棒性较强。Georgescu等人（2002）[55]实验研究了气泡在不同表面张力液体自由表面处的破裂，发现表面张力大的液体形成的射流水滴体积大于表面张力小的液体，且原始气泡对应的临界半径（可形成射流水滴的最大气泡半径）也越大。

在国内，黄继汤等人（1996）[80]研究了表面张力对于单个气泡运动特性的影响，研究表明表面张力在气泡膨胀阶段起到延缓作用，在收缩阶段则起到了加速作用。并通过高速摄影实验对五种表面张力不同的液体内气泡运动进行拍摄分析，分析得到表面张力越小的液体中，气泡壁面速度达到的最大值越大。王起棣和张慧生等人（2003，2005）[81-82]分别研究了表面张力对于刚性壁面附近单个和两个气泡运动的影响，数值模拟采用轴对称边界元法，结果表明表面张力系数越大的流体中气泡溃灭速度越快，射流随之变宽变短。而对于两个气泡情况而言，表面张力的影响与壁面和上气泡对于下气泡的Bjerknes力的大小有关。

4.3 粘性

关于气泡粘性效应的考虑，最直接的做法是求解NS方程，然而由于数值直接求解NS方程的难度，学者们通常从高雷诺数和低雷诺数两方面入手简化问题。对于高雷诺数的情况，最经典的当属Miksis等人（1982）[83]引入的气泡边界层理论，Miksis等人认为气泡周围存在一薄边界层，边界层外流体满足无粘不可压缩势流理论，边界层内考虑流体粘性效应，边界层两侧满足法向应力连续，但不满足切向应力连续。由于此边界层很薄，忽略边界层内外速度的变化，这样便可采用边界积分法数值模拟计及弱粘性效应的气泡动态变化。此后，研究人员从不同角度改进Miksis等人的气泡边界层理论。Lundgren和Mansour（1988）[84]进一步分析边界层方程，获得了边界层两侧压力差和法向速度差的表达式。Boulton-Stone和Blake（1993，1995）[61-62]又进一步扩展了Lundgren和Mansour的方法，使得边界层两侧速度的切向分量也连续，在此基础上首次成功地数值模拟了单个气泡在自由面处的破裂现象。Georgescu等人（2002）[55]也采用Miksis的气泡边界层理论模拟气泡在自由面处的破裂现象，在边界积分和时间步进上均采用二阶单元，以获得更好的精度。Klaseboer等人（2011）[85]采用Joseph和Wang（2004）[86]的思想修正边界层理论，同时考虑了法向和切向应力的连续，应用边界积分法数值模拟了高雷诺数下气泡在流体中的脉动上浮现象。

对于低雷诺数问题，当雷诺数很低时，流场可简化为Stokes流。Youngren和Acrivos（1976）[87]首先应用边界积分方程求解了Stokes流内单个对称气泡的变形，对于表面张力较大的气泡形状，数值计算结果与Barthes-Biesel和Acrivosde（1973）[88]的解析解吻合良好；对于表面张力较小的气泡形状，与Buckmaster（1972）[89]应用细长体理论获得的表达式基本一致。 近期，Eshpuniyani等人（2008，2010）[90-91]采用边界元法研究了Stokes流中贴附于壁面上的微气泡的变形和滑移，旨在为强粘性微血管循环流中气泡的运动形态提供参考。

目前在国内，关于考虑粘性的气泡运动，如粘性流中气泡上升问题等，多采用直接数值模拟法如Front Tracking（陈斌等人，2005[92]），Level Set（李彦鹏等人，2007[93]）等，很少有对边界积分法进行粘性修正的。关于边界积分法考虑粘性问题还有待于进一步研究。

5 需进一步研究的问题

国内外关于气泡的不稳定界面现象的研究日益丰富，实验技术和数值模拟手段相应地完善成熟。尤其对于边界元法，在处理气泡这类具有大变形气液交界面的问题上具有较大的优势，获得了研究人员的青睐，取得了较多有益的进展与结论。不过也存在以下几方面不足：

（1）目前建立的气泡边界元模型中，通常忽略了粘性效应的影响，而在气泡减阻技术或气泡在漩涡流场中的运动等情况，粘性效应可能会很大地影响气泡的运动，并改变气泡的载荷形式，关于如何在边界元模型中考虑粘性修正，仍需要加强研究；

（2）目前关于气泡的撕裂和融合现象主要集中在轴对称模型，关于三维气泡的撕裂和融合数值模型和准则的研究相对薄弱，尚需加强。对于融合问题，同向融合研究成果较多，异相融合的相关数值模拟尚需深入；

（3）气泡和自由面的相互作用研究比较丰富，但关于水下爆炸气泡在接近自由面附近的运动，包括气泡在自由面的破裂和自由面自身的飞溅等研究尚需加强；

（4）目前业内基本已公认气泡射流是造成附近结构毁伤的主要载荷之一，然而气泡射流对结构的毁伤模式，以及不同爆距等参数如何影响气泡射流对结构的毁伤等问题，仍不十分明朗，需要进一步深入研究。

[1]Hsiao C T,Pauley L L.Numerical computation of the tip vortex flow generated by a marine propeller[J].Journal of Fluids Engineering,1999,121(3):638-645.

[2]van Wijingaarden L.On the collective collapse of a large number of gas bubbles in water[C]//Proceeding of 11th International Congress on Applied Mechanics.Munich,1964:854-861.

[3]Li Y B,Wu X Y,Ma Y,Wang J G.A method based on potential theory for calculating air cavity formation of an air cavity resistance reduction ship[J].Journal of Marine Science Application,2008,7:98-101.

[4]Unverdi S O,Tryggvason G.A front-tracking method for viscous,incompressible,multi-fluid flows[J].Journal of Computational Physics,1992,100:25-37.

[5]Bonometti T,Magnaudet J.An interface-capturing method for incompressible two-phase flows.Validation and application to bubble dynamics[J].International Journal of Multiphase Flow,2007,33:109-133.

[6]Hirt C W,Nichols B D.Volume of Fluid(VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics,1981,39:201-225.

[7]Huang J T,Zhang H S.Level set method for numerical simulation of a cavitation bubble,its growth,collapse and rebound near a rigid wall[J].Acta Mechanica Sinica,2007,23:645-653.

[8]Ding Z,Gracewsi S M.The behaviour of a gas caviyy impacted by a weak or strong shock wave[J].Journal of Fluid Mechanics,1996,309:183-209.

[9]Longuet-Higgins M S,Cokelet E D.The deformation of steep surface waves on water.I.A numerical method of computation[J].Philosophical Transactions of the Royal Society of London,1976,A 350:1-26.

[10]Lenoir M.A calculation of the parameters of the high-speed jet formed in the collapse of a bubble[J].Journal of Applied Mechanical Technology Physics,1979,20(3):333-337.

[11]Guerri L,Lucca G,Prosperetti A.A numerical method for the dynamics of nonspherical cavitation bubbles[C]//In Proceedings of the Second International Colloquium on Drops and Bubbles.California,USA,1981:175.

[12]Blake J R,Taib B B,Doherty G.Transient cavities near boundaries.Part I:Rigid boundary[J].Journal of Fluid Mechanics,1986,170:479-497.

[13]Blake J R,Taib B B,Doherty G.Transient cavities near boundaries.Part 2:Free surface[J].Journal of Fluid Mechanics,1987,181:197-212.

[14]Chahine G L,Numerical modeling of the dynamic behavior of bubble in nonuniform flow field[C]//In ASME symposium on numerical methods for multiphase flows.Toronto,1990.

[15]Best J P.The formulation of toroidal bubbles upon collapse of transient cavities[J].Journal of Fluid Mechanics,1993,251:79-107.

[16]Zhang S,Duncan J H.On the nonspherical collapse and rebound of a cavitation bubble[J].Physics of Fluids,1994,6(7):2352-2362.

[17]Brujan E A,Keen G S,Vogel A,et al.The final stage of the collapse of a cavitation bubble close to a rigid boundary[J].Physical of Fluids,2002,14(1):85.

[18]Lee M,Klaseboer E,Khoo B C.On the boundary integral method for the rebounding bubble[J].Journal of Fluid Mechanics,2007,570:407-429.

[19]Bjerknes.Fields of force[M].Columbia:Columbia University Press,1966:45-47.

[20]Blake J R,Gibson D C.Cavitation bubbles near boundaries[J].Journal of Fluid Mechanics,1987,19:99-123.

[21]Wang Q X,Yeo K S,Khoo B C,Lam K Y.Vortex ring modelling of toroidal bubbles[J].Theoretical Computational Fluid Dynamics,2005,00:1-15.

[22]张阿漫,姚熊亮.基于边界积分法的气泡动态特性综述[J].力学进展,2008,38(5):561-570.

[23]张阿漫,汪 玉,闻学友,倪宝玉,等.水下爆炸气泡动态特性研究综述[J].船舶力学,2009,13(5):828-840.

Zhang Aman,Wang Yu,Wen Xueyou,Ni Baoyu,et al.Review of the dynamics of the underwater explosion bubble[J].Journal of Ship Mechanics,2009,13(5):828-840.

[24]Choi J H.Dynamics and noise emission of vortex cavitation bubbles[R].University of Michigan,USA,2006.

[25]Chahine G L,Genoux Ph F.Collapse of a cavitating vortex ring[J].Trans.ASME,Journal of Fluids Engineering,1983,105:400.

[26]Duraiswami R,Chahine G L.Analytical study of a gas bubble in the flow field of a line vortex[J].Cavitation and Multiflow Forum,ASME,1992.

[27]Hsiao C T,Pauley L L.Study of tip vortex cavitation inception using Navier-Stokes computation and bubble dynamics model[J].Journal of Fluids Engineering,1999,l21:198-204.

[28]Hsiao C T,Chahine G L.Scaling effect on prediction of cavitation inception in a line vortex flow[J].Journal of Fluid Engineering,2003,125:53-60.

[29]Hsiao C T,Chahine G L.Prediction of tip vortex cavitation inception using coupled spherical and nonspherical bubble models and Navier-Stokes computations[J].Journal of Marine Science Technology,2004,8:99-108.

[30]Hsiao C T,Jain A,Chahine G L.Effect of gas diffusion on bubble entrainment and dynamics around a propeller[C]//26th Symposium on Naval Hydrodynamics,September,17-22,2006.Rome,Italy,2006.

[31]Oweis G F,Hout I E,Iyer C,Tryggvason G,Ceccioa S L.Capture and inception of bubbles near line vortices[J].Physics of Fluids,2005,17:1-14.

[32]Ni B Y,Zhang A M,Yao X L,Wang B.Numerical simulation of trajectory and deformation of bubble in tip vortex[J].Applied Mathematics and Mechanics(English edition),2012,33(6):1-16.

[33]王献孚.空化泡和超空化泡流动理论及应用[M].北京:国防工业出版社,2009.

[34]韩宝玉,熊 鹰,叶金铭.梢涡空泡初生预报方法[J].海军工程大学学报,2011,23(4):27-32.

[35]Bourne N K,Field J E.Shock-induced collapse of single cavities in liquids[J].Journal of Fluid Mechanics,1992,244:225-240.

[36]Blake J R,Keen G S,Tong R P,Wilson M.Acoustic cavitation:the fluid dynamics of non-spherical bubbles[J].Philosophy Transaction of Royal Society,London,1999,A 357:251-267.

[37]Klaseboer E,Turangan C K,Fong S W,et al.Simulations of pressure pulse-bubble interaction using boundary element method[J].Computer Methods in Applied Mechanics and Engineering,2006,195:4287-4302.

[38]Klaseboer E,Fong S W,Turangan C,et al.Interaction of lithotripter shockwaves with single inertial cavitation bubbles[J].Journal of Fluid Mechanics,2007,593:33-56.

[39]倪宝玉,张阿漫,朱 枫,吴巧瑞.舰船近场爆炸反射波对气泡射流特性影响[J].哈尔滨工程大学学报,2011,32(8):976-983.

[40]Jamaluddin A R,Ball G J,Leighton T G.Free-lagrange simulations of shock/bubble interaction in shock wave lithotripsy[C].Proceedings of the Second International Conference on Computational Fluid Dynamics,2002.

[41]Kim J O,Kim S D.Bubble breakage phenomena,phase holdups and mass transfer in three-phase fluidized beds with floating bubble breakers[J].Chemical Engineering Process,1990,28:101-111.

[42]Wilkinson P M,Schayk A V,Spronken J P M.The influence of gas density and liquid properties on bubble breakup[J].Chemical Engineering Science,1993,48(7):1213-1226.

[43]Martinez-Bazan C,Montanes J L,Lasheras J C.On the breakup of an air bubble injected into a fully developed turbulent flow.Part 1.Breakup frequency[J].Journal of Fluid Mechanics,1999,401:157-182.

[44]Ishida H,Nuntadusit C,Kimoto H.Cavitation bubble behavior near solid boundaries[J].CAV,2001,A5-003.

[45]Choi J K,Chahine G L.Noise due to extreme bubble deformation near inception of tip vortex cavitation[J].FEDSM,2003:45313.

[46]Ellis A T.Techniques for pressure pulse measurements and high-speed photography in ultrasonic cavitation,cavitation in hydrodynamics[C].HMSO,London,Paper 8,1956.

[47]Lauterborn W,Hentschel W.Cavitation bubble dynamics studied by high speed photography and holography.Part I:Ultrasonics[M].In press,1985:258-293.

[48]Mohammad T S,Amir A,Mohammad R S.Numerical study on the dynamics of an electrical discharge generated bubble in EDM[J].Engineering Analysis with Boundary Elements,2006,30:503-514.

[49]张阿漫,倪宝玉,宋炳月,姚熊亮.狭窄流域内气泡破裂现象数值模拟[J].应用数学和力学,2010,31(4):420-432.

[50]Rungsiyaphornrat S,Klaseboer E,Khoo B C,Yeo K S.The merging of two gaseous bubbles with an application to underwater explosions[J].Computers&Fluids,2003,32:1049-1074.

[51]Chen P,Sanyal J,Dudukovi M P.Numerical simulation of bubble columns flows:effect of different breakup and coalescence closures[J].Chemical Engineering Science,2005,60:1085-1101.

[52]Chevaillier J P,Klaseboer E,Masbernat O,Gourdon C.Effect of mass transfer on the film drainage between colliding drops[J].Colloid and Interface Science,2006,299:472-485.

[53]Klaseboer E,Chevaillier J P,Gourdon C,Masbernat O.Film drainage between colliding drops at constant approach velocity:Experiments and modeling[J].Colloid and Interface Science,2000,229:274-285.

[54]张阿漫,曾令玉,王诗平,杨树涛,姚熊亮.水下爆炸气泡融合动态特性研究[J].应用数学和力学,2010:31(2):163-170.

[55]Georgescu S C,Achard J L,Canot E.Jet drops ejection in bursting gas bubble processes[J].European Journal of Mechanics B/Fluids,2002,21:265-280.

[56]Malysa K,Krasowska,Krzan M.Influence of surface active substances on bubble motion and collision with various interfaces Advances in Colloid and Interface Science,2005,114-115:205-225.

[57]Kientzler C F,Arons A B,Blanchard D C,Woodcock A H.Photographic investigation of the projection of droplets by bubbles bursting at a water surface[J].Tellus,1954,6(1):1-7.

[58]Garner F H,Ellis S R M,Lacey J A.The size distribution and entrainment of droplets[J].Transaction of Institute Chemistry Engineering,1954,32:222-235.

[59]Resch F J,Afeti G M.Gas-liquid conversion at a free interface and its environmental implications[J].Fluid Mechanics and Thermodynamics,1991:1066-1076.

[60]Afeti G M,Resch F J.From breaking waves to breaking bubbles:Production characteristics of the marine droplet aerosol[J].Trends Geophysics Research,1992,1:113-125.

[61]Boulton-Stone J M,Blake J R.Gas bubbles bursting at a free surface[J].Journal of Fluid Mechanics,1993,254:437-466.

[62]Boulton-Stone J M.The effect of surfactant on bursting gas bubbles[J].Journal of Fluid Mechanics,1995,302:231-257.

[63]Sussman M,Smereka P.Axisymmetric free boundary problems[J].Journal of Fluid Mechanics,1997,341:269-294.

[64]Prosperetti A,Lezzi A.Bubble dynamics in a compressible liquid.Part I.First-order theory[J].Journal of Fluid Mechanics,1986,168:457-478.

[65]Herring C.Theory of the pulsations of the gas bubble produced by an underwater explosion[R].Office of Science Research and Development Report 236(NDRC Report C-4-sr 10-010,Columbia University),1941.

[66]Keller J B,Kolodner I I.Damping of underwater explosion bubble oscillations[J].Journal of Applied Physics,1956,27:1152-1161.

[67]Lezzi A,Prosperetti A.Bubble dynamics in a compressible liquid.Part II.Second-order Theory[J].Journal of Fluid Mechanics,1987,185:289-321.

[68]Pozrikidis C.Expansion of a compressible gas bubble in Stokes flow[J].Journal of Fluid Mechanics,2001,442:171-189.

[69]Darren G,Crowdy.Compressible bubbles in Stokes flow[J].Journal of Fluid Mechanics,2003,476:345-356.

[70]Wang Q X,Blake J R.Non-spherical bubble dynamics in a compressible liquid.Part 1:Travelling acoustic wave[J].Journal of Fluid Mechanics,2010,659:191-224.

[71]Wang Q X,Blake J R.Non-spherical bubble dynamics in a compressible liquid.Part 2:Acoustic standing wave[J].Journal of Fluid Mechanics,2011,679:559-581.

[72]Geers T L,Hunter K S.An integrated wave-effects model for an underwater explosion bubble[J].Acoustic Society of America,2002,111(4):1584-1601.

[73]王诗平.水中结构物附近爆炸气泡运动特性研究[J].哈尔滨:哈尔滨工程大学,2011.

[74]Beeching R.Resistance to cavitation erosion[J].Transaction Institute of Engineering and Shipbuilders,1942,85:273-277.

[75]Akhatov I,Gumerov N,Ohl C D,Parlitz U,Lauterborn W.The role of surface tension in stable single-bubble sonoluminescence[J].Physical Review Letters,1997,78(2):227-230.

[76]Bronshtein I N,Semendyayev K A.Handbook of Mathematics[M].3rd edition.Springer Publication,Berlin,1997.

[77]Brackbill J U,Kothe D B,Zemach C.A continuum method for modeling surface tension[J].Journal of Computational Physics,1992,100:335-354.

[78]Zinchenko A Z,Rother A M,Davis R H,A novel boundary-integral algorithm for viscous interaction of deformable drops[J].Physics Fluids,1997,9(6):1493-1511.

[79]Rallison J M.A numerical study of the deformation and burst of a viscous drop in general shear flows[J].Journal of Fluid Mechanics,1981,109:465.

[80]黄继汤,陈嘉范,丁 彤,田立言.表面张力对单空泡运动特性的影响[J].水利学报,1996,12:1-7.

[81]王起棣,张慧生.表面张力对固壁近旁空化气泡溃灭特性的影响[J].复旦学报（自然科学版）,2003,42(2):201-207.

[82]张振宇,王起棣,张慧生.表面张力对近固壁二空化泡影响的数值研究[J].力学学报,2005,37(1):100-104.

[83]Miksis M J,Vanden-Broeck J M,Keller J B.Rising bubbles[J].Journal of Fluid Mechanics,1982,123:31-41.

[84]Lundgren S,Mansour N.Oscillations of drops in zero gravity with weak viscous effects[J].Journal of Fluid Mechanics,1988,194:479-510.

[85]Klaseboer E,Manica R,Chan D Y C,Khoo B C.BEM simulations of potential flow with viscous effects as applied to a rising bubble[J].Engineering Analysis with Boundary Elements,2011,35:489-494.

[86]Joseph D D,Wang J.The dissipation approximation and viscous potential flow[J].Journal of Fluid Mechanics,2004,505:365-377.

[87]Youngren G K,Acrivos A.On the shape of a gas bubble in a viscous extensional flow[J].Journal of Fluid Mechanics,1976,76(3):433-442.

[88]Barthes-Biesel D,Acrivosde A.Deformation and burst of a liquid droplet freely suspended in a linear shear field[J].Journal of Fluid Mechanics,1973,61:1.

[89]Buckmaster D.Pointed bubbles in slow viscous flow[J].Journal of Fluid Mechanics,1973,55:385.

[90]Eshpuniyani B,Fowlkes J B,Bull J L.A boundary element model of microbubble sticking and sliding in the microcirculation[J].International Journal of Heat and Mass Transfer,2008,51:5700-5711.

[91]Calderon A J,Eshpuniyani B,Fowlkes J B,Bull J L.A boundary element model of the transport of a semi-infinite bubble through a microvessel bifurcation[J].Physics of Fluids,2010,22,061092(1-11).

[92]陈 斌,Kawamura T,Kodama Y.静止水中单个上升气泡的直接数值模拟[J].工程热物理学报,2005,26(6):980-982.

[93]李彦鹏,关卫省.应用Level Set方法模拟三维气泡与液滴的运动[J].地球科学与环境学报,2007,29(2):217-220.

Review on the numerical simulation methods of the instability of bubble interface under complex flow field

NI Bao-yu,ZHANG A-man,ZHANG Zhong-yu,LI Shuai

(School of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China)

The study on the instability of the bubble interface is a hot topic.Among the various numerical methods,Boundary element method(BEM)is widely adopted to simulate the bubble dynamics due to its high efficiency.The development and application of BEM and relative numerical treatments on bubble dynamics are reviewed in this paper.The collapse of bubble near various boundaries under various ambient flows is summarized;The splitting,mergence and breakup of the bubble near free surface are concluded and discussed;the influences of three significant factors,namely compressibility,surface tension and viscosity,are reviewed and analyzed;and,the problems needed to be studied further are put forward.

bubble;interface instability;BEM;split;mergence

O35

A

10.3969/j.issn.1007-7294.2014.04.015

1007-7294（2014）04-0470-14

2013-08-03

国家自然科学基金重点项目（50939002）；国家自然科学基金委员会—中国工程物理研究院联合基金资助项目（10976008）；国防基础科研计划（B2420110011）；国家自然科学基金青年项目（51009035）；黑龙江省自然科学基金（E201047,A200901）资助课题；中央高校基本科研业务费专向资金资助（HEUCF140116）

倪宝玉（1986-），男，博士，哈尔滨工程大学船舶工程学院讲师，E-mail:baoyuni@gmail.com；

张阿漫（1981-），男，教授，博士生导师。