Chinese Journal of Ship Research
近冰面航行潜艇阻力及绕流场分析
引用格式：柏铁朝,许建, 王国栋, 等.近冰面航行潜艇阻力及绕流场分析[J]. 中国舰船研究, 2021, 16(2): 36–48. BAI T C, XU J, WANG G D, et al. Analysis of resistance and flow field of submarine sailing near the ice surface[J]. Chinese Journal of Ship Research, 2021, 16(2): 36–48.近冰面航行潜艇阻力及绕流场分析
柏铁朝*1，许建1，王国栋1，余琨1，胡笑寒2 1中国舰船研究设计中心，湖北武汉 430064 2哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 150001
摘 要:［目的］采用数值模拟方法探究潜艇近冰面航行时的水动力性能。［方法］选取Suboff 全附体潜艇模型为研究对象，在STAR-CCM+软件中采用 RANS 方法，再结合SST k-ω湍流模型和体积分数法计算艇体的水动力性能。然后进行网格无关性验证以及近冰面水动力特性分析的方法验证，以确定计算方法的有效性。［结果］结果表明，潜艇近冰面航行时，潜艇的总阻力系数CT在同一弗劳德数下随潜深的增加而减小，在同一潜深下随弗劳德数的增加而减小；在同一弗劳德数下，当无量纲深度大于1.63 时，潜深几乎不影响艇体的阻力系数。［结论］研究表明，潜艇的总阻力系数受弗劳德数以及潜深的影响较大，随着潜深的增大，潜深对阻力系数的影响将逐渐减小。关键词：潜艇；冰面；水动力；扰流场；数值模拟中图分类号: U661.31+1 文献标志码:A DOI：10.19693/j.issn.1673-3185.01975
Analysis of resistance and flow field of submarine sailing near the ice surface
BAI Tiechao*1, XU Jian1, WANG Guodong1, YU Kun1, HU Xiaohan2 1 China Ship Development and Design Center, Wuhan 430064, China 2 College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China
Abstract: ［Objectives］This paper studies the hydrodynamic performance of the submarine sailing near the ice surface using numerical simulation method.［Methods］the Suboff full-appended model, RANS method, SST k-ω turbulence model and volume of field method are used to calculate the hydrodynamic performance of the submarine in the STAR-CCM+ software. The effectiveness of the calculation method is guaranteed by the method of grid independence verification and the method of hydrodynamics analysis near the ice surface.［Results］The results show that when the submarine sails under ice surface, the total resistance coefficient of submarine decreases with the increase of the submergence depth at the same Froude number, and decreases with the increase of Froude number at the same depth. In addition, when the dimensionless depth is greater than 1.63, the submergence hardly affects the drag coefficient.［Conclusions ］ Studies have shown that the total drag coefficient of submarine is greatly influenced by Froude number and submergence depth. With the increase of submergence depth, the influence of submergence depth on drag coefficient decreases gradually.
Key words: submarine；ice surface；hydrodynamics；flow field；numerical simulation
0 引 言
作为海军的重要作战力量，潜艇的隐蔽性使其具有强大的威慑力，其探测与反探测研究一直受到海军强国的高度重视[1]。极地环境的复杂性使得潜艇在冰区的航行情况与自由水域相比差异
较大。冰面的覆盖改变了水面处的边界条件，并且考虑到潜艇航行时产生的挠曲重力波，有可能导致冰面破碎，从而使航行情况更为复杂[2]。因此，研究挠曲重力波对冰面的影响，对执行军事任务具有重要意义，对潜艇近冰面航行的水动力特性分析也十分重要。在潜艇水动力分析研究初期，主要是致力于无界黏性流域条件下的情况。截至20 世纪 80 年代中期，许多国家进行了大量的数值研究，其在原有势流计算方法和切片理论的基础上予以发展，考虑了自由液面带来的影响[3]。至 20 世纪 80年代后期，美国国防先期研究计划局（ DARPA）发起建立了带各种典型附体的轴对称回转体流场试验数据库[4]。近年来，针对潜艇近自由液面航行的研究开始逐渐增多，并以数值研究为主。Jagadeesh 等[5] 基于 RANS 方法，讨论了自由液面对回转体潜艇模型航行的影响。Saout 等[6] 针对潜艇近自由液面航行时的稳定性进行了研究。Dawson[7]就浸没深度、速度和细长比对流线型对称体潜艇与自由液面的相互作用开展了相关试验和数值研究。我国对潜艇水动力性能的研究主要集中在对 Suboff 潜艇全附体模型阻力、升力、表面压力及精细流场的数值计算上。李艳等[8] 基于 RANS方法，将湍流模型与壁面函数法相结合，对Suboff全附体潜艇模型进行了数值模拟，并将计算值与试验值进行对比，验证了模拟方法的准确性。张楠等[9-10] 基于 RANS 方法，结合流体体积（volume of fluid，VOF）法，分别计算了 Suboff 近自由液面航行时的阻力以及各潜深下的增阻系数、艇−桨干扰特性，并讨论了围壳对自由液面兴波的影响。邱云明和邓锐[11] 对 Suboff全附体潜艇模型近水面航行时的兴波进行了多工况下的计算，并对比了各弗劳德数下的自由液面波形和中纵剖面上的波高。然而，无论国内还是国外，少有关于Suboff 潜艇模型在冰面下航行时的流场特性及其规律研究。而对冰区潜艇近冰面航行进行数值模拟，掌握其阻力与流场特性进而探索挠曲重力波破冰方法意义重大。为了系统分析潜艇近冰面航行时的阻力及流场分布规律，本文拟以 Suboff 全附体潜艇模型为研究对象，在 STAR-CCM+软件中应用RANS 方法计算艇体的水动力性能，然后对不同潜深、航速等条件下艇体近冰面航行时的阻力系数、表面压力系数、表面摩擦系数以及艇体绕流场特性进行分析。
1 数值方法1.1 基本控制方程
对于计算流体力学问题，相关的控制方程包括质量守恒方程和动量方程。质量守恒方程：
式中：ρ 为流体密度；ui为速度矢量在3个方向的速度分量；xi为基于速度分量的3个方向的位移； u为速度矢量；p 为压强；Fi 为质量力；t 为时间； τxi，τyi，τzi 为黏性应力τ的分量。
1.2 湍流模型
本文所使用的湍流模型为SST k-ω 模型[12-13]。其中湍流动能k的输运方程为：
式中：xj 为基于湍流动能和耗散率的3 个方向的位移；Gk 为由层流速度梯度产生的湍流动能；Gω为由ω方程产生的湍流动能；Tk和 Tω分别为k和ω的扩散率；Yk 和 Yω 为由扩散产生的湍流；Dω为正交发散项。
2 物理条件设置2.1 计算对象参数
本文以美国DARPA 提出的 Suboff 全附体潜艇模型为研究对象，其几何模型及参数分别如图 1和表1所示。
2.2 计算域及边界条件设置
由于所采用的 Suboff 全附体潜艇模型和流域
为对称的，所以通过设置对称面计算一半的流域即可。潜艇近冰面航行时，在潜艇穿过冰盖下某一位置的瞬间，可以认为冰盖尚未变形，或者是其变形对潜艇绕流场的影响微弱，因此，将艇体上方边界设置为无滑移壁面边界。为了模拟无限流域的物理背景，将来流出口设置为压力出口，潜艇中纵剖面设置为对称面，潜艇表面设置为无滑移壁面。为保证计算稳定性，其余边界设置为速度进口。具体边界条件的设置如图2 所示。图中，L为艇体长度，d为模型潜深大小。
2.3 网格收敛性分析
本文采用3种不同的网格尺寸进行验证。在划分网格时，通过调整网格的基本尺寸，可以在不改变网格相对尺寸的条件下，改变网格的疏密程度。对艇体加密时，选择面网格加密方法，并将艇体表面分割为指挥室围壳、尾翼和艇体3部分进行计算。采用的网格分为3套，其基础尺寸分别为 0.8，1.0 和 1.2 m，各基础尺寸下的网格数量信息（半侧计算域）如表2所示。边界层取10 层， Y+值为30，考察艇体前进方向所受的阻力。计算域的网格划分如图3所示。目前，直航全附体潜艇阻力测试普遍参照的是 1998年美国海军水面战中心卡德洛克分部
（NSWC） [14]提出的不同航速下所对应的总阻力试验值。选择航行速度 V=3.05 和 5.14 m/s 进行网格无关性验证时，各网格尺度下的计算结果及误差如表3和表4所示。无关性仿真计算结果与实体试验的比较说明网格尺寸为单调收敛，将选择误差最小的基础尺寸 0.8 m作为 Suboff 潜艇模型变参数仿真计算的
网格划分参数。
2.4 计算方法验证
为了验证本文数值模拟方法的可靠性，设置了多个航行速度，用以计算无界绕流条件下Suboff全附体潜艇模型的直航阻力值，并与文献[4] 中的试验数据进行了对比，结果如图4所示。由图4可以看出，仿真值和试验值均是随航速的增加而增加的，二者吻合较好，仿真值只在较高航速时略高于试验值，整体来说误差控制在2%以内，验证了本文采用方法的可行性。
2.5 近冰面航行计算模型建立方法验证
由于基于重叠网格的方法（方法1），也即采用潜艇在静止的冰面下直线航行的计算模型开展数值计算的方法，存在计算资源耗费大、计算速度过慢及计算结果不易收敛的问题，所以在实际计算中采用的是潜艇在均匀来流中的水动力性能分析方法（方法2）来模拟艇−冰−流之间的相互作用，由此避免重叠网格的应用，降低计算难度并提高计算速度。这种方法的几何模型及网格划分示意图如图5所示。但是，采用艇不动的方法（方法2）会消除艇与冰面之间的相对运动。针对这一缺陷，分别提取基于方法1和方法2 这 2种方法计算得到的潜艇航行阻力、艇体表面压力分布、中纵剖面上半缘线压力系数分布及冰面下表面压力分布结果进行了对比验证，以此确认方法2中计算模型设置对计算结果的影响程度。各计算结果的对比如表 5 和图6、图 7所示。计算工况选择较高航速V=5.29 m/s（潜艇以较高航速航行时对计算方法的准确性及稳定性要求较高，作为方法1/方法 2的对比验证计算工况较为合理）下的不同潜深（d= 1.67， 2.08 m ）。观察表5，发现在不同潜深下，相比于方法1，方法 2在潜艇阻力计算结果上的偏差均在3%以内。
为了更加精确地对比2种计算方法对艇体表面压力的影响，提取了各方法下艇体中纵剖面上半缘线的表面压力系数CP 分布，如图6 所示（图中，横坐标为艇不同横向位置的无因次值）。由图可见，基于方法1 和方法2得到的艇体中纵剖面表面压力系数分布基本相同，仅在围壳后端、尾翼、艇艉型线收缩处出现了略大的差异。与方法 1的计算结果相比，基于方法2计算的艇体中纵剖面上半缘线表面压力系数偏差在1.6% 以内。图 7 所示为航速V =5.29 m/s、潜深d =1.67 m时，基于2种计算方法得到的冰面下表面压力分布云图。其中，横坐标为艇不同横向位置的无因次值，纵坐标为2种计算方法所得到的压力系数的物理值。由图可见，采用2种计算方法得到的(b) 艇体加密
冰面下表面压力分布形态基本一致，因艇体与冰面之间的水流流速较快，导致冰面在艇体上方出现了环状低压区，并在低压区的前、后方出现了环状高压区，这种压差在一定程度上会导致冰板的形变。对比图7(a) 与图 7(b) 可见，采用方法2模拟潜艇近冰面航行冰面所受到的影响时，计算模型的设计缺陷并未对结果形成较大的干扰。综上，考虑到采用方法1开展潜艇近冰面航行计算资源耗费较大、计算速度过低，因此在有限的计算资源与时间条件下，本文将选择方法2作为潜艇近冰面航行数值模拟的计算方法。对于方法2 在计算模型设置上存在的缺陷，通过对潜艇航行阻力、艇体表面压力系数及冰面下表面压力分布计算结果的对比验证，证明此缺陷对近冰面航行潜艇的阻力、绕流场及冰面下表面压力的数值模拟结果影响较小，不影响方法2 的可行性。
2.6 计算工况设置
参考 Pogorelova 等[15] 对潜艇在弹性板下运动的试验研究，数值模拟的工况在实尺度上与该试验中的工况设置一致，如表6 所示。表中，Fr为艇长弗劳德数，d为 Suboff 潜艇航行深度（艇体中轴线至冰面下表面的距离），D为艇体最大直径（本文计算模型取 D=0.508 m）。
3 计算结果与分析3.1 阻力预报结果分析
各工况下近冰面航行潜艇总阻力（单位：N）计算结果如表7所示，各潜深下艇体的各成分阻力系数随弗劳德数Fr的变化规律如图8所示。其中，摩擦阻力系数CF，压阻力系数 CVP 和总阻力系数CT的定义如下：
式中：RF，RPR，RT 分别为艇体的摩擦阻力、压阻力和总阻力；ρ=998.16 kg/m3，为水密度；A 为艇体湿表面积，由于本数值模拟时潜艇是完全浸入水中的，所以A为 Suboff 全附体潜艇模型的表面积， A=6.361 m3。由图 8(a) 可知，潜艇的总阻力系数CT 在同一弗劳德数Fr 下随潜深d的增加而降低，在同一潜深 d下随弗劳德数Fr的增加而降低。并且，在不同潜深下，艇体阻力系数随航速的变化趋势几乎相同。潜艇近冰面航行时，总阻力包括摩擦阻力和压阻力。对比图8(b) 和图 8(c) 中的数值发现，潜艇近冰面航行时，在设计工况范围内，摩擦阻力系数 CF占总阻力系数CT的主要部分。通过对潜艇周围阻力性质的分析不难得出，由于冰面的存在，使得艇体周围不存在兴波影响，因而减小了艇体所受到的压阻力；此外，由于艇体全部浸没于水中，且潜艇的湿表面积较大，故进一步削减了压阻力对总阻力的贡献。摩擦阻力为边界层内剪切应力的合力。由于壁面剪切应力τ在湍流流态时正比于航速V 3/2 ，所以壁面剪切应力τ随V的增大情况小于V 2关系，从而导致了局部摩擦阻力系数Cτ 随航速的升高而降低的情况[16]。因为 Cτ 与 CF的变化规律是一致的，所以，CF随航速的升高有降低趋势，这与图8(c) 中 CF 随 Fr 的变化趋势是一致的。在各个潜深下，4条曲线几乎重合，这说明在潜艇近冰面航行情况下，潜深对CF的影响很小，这一结论在文献 [15]中得到了印证。观察图 8(b) 并与摩擦阻力系数CF 相比，发现压阻力系数CVP随潜深的变化较大。随着艇体逐渐远离冰面，CVP 将逐渐降低，但随着潜深的增大，压阻力系数CVP的变化趋势逐渐减小。具体体现为，当潜深 d=1.25 ，1.67 和 2.08 m（d/D=2.46， 3.29，4.09 ）时，所得到的艇体CVP 值几乎没有变化，不过在最小潜深 d=0.83 m（d/D=1.63）时 CVP 的变化相对较大。所以，在近冰面航行时，当 d/D>1.63