潜艇阻力及流场数值仿真策略优化分析

李 鹏， 王 超， 孙华伟， 郭春雨

(哈尔滨工程大学 船舶工程学院， 哈尔滨 150001)

随着计算机的发展，计算流体力学(Computational Fluid Dynamics, CFD)被越来越多地应用到海洋结构物的性能研究中.潜艇是重要的海洋结构物，准确了解其不同航速下的水动力性能、精确掌握其流场特征，对提高潜艇自身性能具有重要意义[1].

国内外学者针对潜艇性能开展过大量的数值研究.Pan等[2]采用RANSk-ω(k为湍动能，ω为单位湍动能耗散率)和雷诺应力湍流模型(Reynolds Stress Model, RSM)对潜艇的阻力进行预报，取得了较好的结果.Posa等[3-5]基于大涡模拟湍流模型(Large Eddy Simulation, LES)对潜艇的尾流场进行仿真计算，结果表明LES能很好地捕捉流动特征，同时为流噪声的预报提供精确输入.Fu等[6-8]就分离涡湍流模型(Detach Eddy Simulation, DES)下的潜艇水动力性能及流场特征进行仿真，结果证实前者对垂直壁面方向的无量纲化距离y+的敏感性微弱，后者对空间离散精度的高需求性.此外，李士强等[8]还就k-ω、k-ε(ε为湍流耗散率)模型对潜艇附体变动的影响进行仿真，并证实关于舵失速的问题，在失速前k-ε获取的涡流脱落频率较k-ω更高.

综上，国内外学者就直航潜艇的性能进行较多分析，但未阐明不同湍流模型在潜艇水动力性能及流场特征计算中的优势或短板.本文针对这一问题，首先展多套网格及时间步长的收敛性分析.接着，针对不同的湍流模型开展特定工况的数值仿真计算，比较预报结果的差异.最后，基于LES-Smagorinsky模型的预报结果对潜艇涡系特征进行分析.

1 理论基础

1.1 控制方程

时均化后不可压缩流体的连续性方程与动量方程为

(1)

(2)

1.2 湍流模型

由于在时间和空间上对计算成本的需求较高，目前主要的CFD求解器均采用求解平均量或滤波量的方法来模拟小尺度流动以节省计算成本.本文就多个湍流模型开展潜艇的水动力性能及流场特征计算，对比各种湍流模型的计算结果，进而筛选适合潜艇水动力性能及流场特征仿真计算的湍流模型.

雷诺平均纳维尔-斯托克斯(Reynolds-Averaged Navier-Stokes, RANS)是一类湍流模型集合的总称，涉及多个子模型，这些模型的特征都是在N-S(Navier-Stokes)方程的基础上添加了一个额外的黏度项且都包含一个通用项湍动能k，即单位质量湍流脉动的动能，而这些模型的区别在于对额外的黏度项处理方式不同：SA(Spalart-Allmaras)模型是一方程湍流模型，无壁面函数要求但增加一个新的变量SA黏度，其优点是对计算内存要求较低，有良好的收敛性，缺点是不适合求解剪切流与分离流；k-ε模型基于k和ε建立，是具有良好收敛性的两方程模型，自提出后逐渐演化出realizablek-ε，RNGk-ε等，局限性体现在对无滑移壁面、逆压梯度等的流动模拟不精确，耗散不易计算；k-ω模型计算结果对初始条件较为敏感且收敛相对困难，但相对于ε，ω更容易求解.LES是基于自相似理论，通过过滤函数将大涡和小涡分离开，大涡直接解析，小涡是各项同性，因此可对其进行模化计算，在解析大涡和模化小涡的过程中，产生的亚格子尺度应力项通过亚格子(Subgrid Scale, SGS) 模型计算.DES模型结合RANS/LES模型各自的优点，在近壁面采用RANS进行计算而利用LES模型的算法对湍流核心区域计算，经过多年发展逐渐演变出延迟分离涡湍流 (Delayed Detach Eddy Simulation, DDES)、改进的延迟分离涡湍流 (Improved Delayed Detached Eddy Simulation, IDDES) 等模型变体，对非定常湍流流动的模拟有较大的优势.

本文基于植入到STAR-CCM+中的湍流模型：RANSk-ω、k-ε、SA、LES-Samgorinsky(为行文简洁，后续以LES-S替代)、WALE(LES-W)、Dynamic Smagorinsky(LES-DS)以及DES-SA、DES-EBKE、DDES、IDDES模型开展Sub-off水动力性能及流场特征的计算，进而甄选合适的湍流模型开展潜艇水动力性能及流场特征的计算.

2 仿真基础

2.1 坐标系和计算域

本文采用的几何模型是Sub-off模型，详细的模型资料可参考文献[9-11]，后续的研究中，在原始几何模型尺度的基础上进行缩放，缩放比λ=0.688 7，模型总长L=3 m，最大直径D=0.35 m，若无特殊说明本文的计算结果均是基于这一参数下的几何模型.

本文设定笛卡尔坐标系Oxyz，其原点位于艇首，如图1(a)所示，图中Uref为自由流场的自由来流速度.本文设定计算域为圆柱，具体尺寸及边界条件的设置见图1(b)，参考文献[12]的研究，当目标几何模型在计算域流向投影面积小于1%时，数值计算结果可忽略阻塞效应的影响，通过简单的数据计算可知，Sub-off模型在计算域上下底面的投影面积占比为0.7%(小于1%)，故可忽略阻塞效应的影响.

图1 几何模型及计算域

2.2 网格划分和工况设定

图2 网格方案3简图

为统一和方便作者进行数值和试验研究，结合计算和实际条件需求，对航速为2.53 m/s的潜艇模型的水动力性能及流场特征进行收敛性计算，以对比不同湍流模型获得的仿真结果.

3 收敛性分析

表2 4套网格方案的计算结果

网格收敛率K=εi, i+1/εi, i-1,εi, i+1是较密两方案结果差，εi, i-1是较稀两方案结果差.K的计算结果可能是K>1，0

表3 网格收敛性计算结果

P=ln (εi, i+1/εi, i-1)lnδ

评估3套网格的收敛指数：

Ii, i+1=DS|ei, i+1|/(δP-1)

式中：DS=1.25为网格收敛安全系数[18]；ei, i+1=εi, i+1.理论上，当P=2时，网格正交性或湍流模型非线性对计算结果准确性的影响程度最小[14,18]，选用中等密度网格进行数值计算具有较高可靠性.可知方案2～4具有更好的收敛率且离散精度更接近理论值，故后续的计算中采用网格方案3进行计算.

表4是不同时间步长计算得到的判定系数结果，以0.001 s的计算结果为基准进行对比.可以看出, 不同时间步长对计算结果的影响较小，为提高计算效率，后续的计算中均采用0.001 s为时间步长.

表4 不同时间步长计算结果

4 结果与讨论

4.1 水动力性能

不同湍流模型对近壁面流动、黏性应力等的解析不同，造成仿真得到的阻力有所差异，表5为不同湍流模型计算得到的潜艇阻力.

表5 不同湍流模型的艇体阻力仿真结果

表6 不同湍流模型计算结果与实验结果对比

4.2 流场特征

采用Q标准来判定流体为涡旋或应变主导[19]，图3所示为不同湍流模型仿真计算得到的Sub-off周围的涡场.图3清晰地表明LES-W、LES-DS对潜艇重要涡特征——马蹄涡的捕捉不足，不同的是除此外的8个湍流模型均很好地捕捉到了围壳附近马蹄涡的流动特征，另外这两种湍流模型预报的潜艇尾肩部的涡结构亦没有其他模型强烈.选用的10种湍流模型中，RANS类湍流模型对潜艇尾部涡强度的预报结果显然低于其他模型，这可能是RANS速度平均后的结果，是RANS模型的算法特点所决定.值得注意的是，10个模型均没能很好地捕捉或预报围壳稍涡的演化，但对于尾舵稍涡的演化均能较好地捕捉.推测原因：一是围壳稍涡的强度没有尾舵稍涡强度大，二是空间分辨率或时间分辨率不足以捕捉其演化过程.限于计算资源以及本文研究重点，在此不做深入探讨.

图3 不同湍流模型计算艇体周围Q采样体积(Q=5 s-2)

伴流是考核潜艇综合性能的重要指标，对伴流的准确预报十分重要.图4所示为不同湍流模型计算得到的Sub-off尾流伴流分数δ随距离变化的曲线，图中，0为潜艇中轴线的最尾端，D、2D、3D、4D为沿着x正方向与潜艇中轴线最尾端的距离.可见，LES-DS、LES-W的预报结果与其他模型有明显偏差，但同属大涡模拟类中LES-S的预报结果有良好的“随群”性，结合对于宏观力的预报结果可知，LES-S模型在预报潜艇宏观变量方面有着较好的适用性.

图4 不同湍流模型预报尾流伴流(0位于艇体尾端)

为了深入探究不同湍流模型对Sub-off流场的预报结果， 本文提取了如图5虚线所示位置的流场时均变量并加以分析，时均流场的采样总长是流经一倍艇长时间. 图6所示为时均轴向速度的周向变化曲线，图中θ为图5中提取截面沿着白色虚线环向位置的角度，u为流场局部轴向速度，通过u/Uref对轴向速度进行无因次化对比, 分析其变化规律，为了有效说明流场，将Sub-off周围流场分为内流场(r/D<0.501)和外流场(r/D>0.501)，其中r为对应曲线的半径值.可知，10种湍流模型对Sub-off的外流场及其平行中体内流场的预报结果吻合程度非常高，推测由于远离壁面或壁面流动简单，各湍流模型对其处理无明显差异，但即使对平行中体内流场的流动预报差异微弱，Sub-off尾内流场的预报差异却十分明显.这表明，湍流模型对因为几何特征变化引发流动变化更为敏感，对流场“传递过程”引发的变化敏感性较弱.图7所示为流场时均压力的周向分布曲线，可以看出，Sub-off周围流场的时均压力系数Cp(Cp=p/(0.5ρUref))分布曲线对湍流模型地变化不敏感，各个模型均有良好的“随群”性，由于尾部附体以及尾部几何收缩造成的负压力梯度，在x/L=0.9处，流场压力波动仍较为剧烈.图6和7的结果表明：时均速度较时均压力对湍流模型的改变更敏感，尤其是近壁面流动，如图6(a)、7(a)所示，这是不同湍流模型调用不同壁面解析函数造成的.图6(b)、7(b)中不同湍流模型得到的流场流速和压力的周向分布基本相同，说明在远流场各个湍流模型预报的数值仿真结果基本相同，但图7(b)中x/L=0.90处的流速和压力分布仍有相对明显的波动，这是由于艇体尾部逆压力梯度造成尾部流速降低，边界层增厚，从而改变流场流速和压力基本对应的变化趋势.

图5 流场提取截面

图6 流场时均轴向速度的周向分布

图7 Cp的周向分布

据上推断，10种湍流模型在Sub-off流场时均变量的捕捉上存在差异，且主要应该体现在流场黏性上，进而导致其他流场变量的变化，在本文采用的时间及空间分辨率的情况下，LES-DE、LES-W同比效果更差，原因可能是空间或时间分辨率不足导致的湍流模型对流场黏性特征的预报存在差异.

不同湍流模型对潜艇周围流场参数的平均量预测结果存在差异，同时对流场的二次量如雷诺应力的预测结果亦有不同.图8所示为不同湍流模型仿真得到的潜艇周围流场的雷诺正应力RNS，其余湍流模型的结果由于数值过小，出于文章简洁的考虑未在文中列出.如同文献[17]，LES模型对流场二次特征的捕捉较好，DES模型次之，RANS模型最差，故研究流场二次量时模型的选择会引起数值仿真结果的巨大差异.

综上可知，就模型尺度下的潜艇而言，RANS类湍流模型不能有效地捕捉潜艇流场的流动细节，对湍流应力有明显的抑制；DES类湍流模型普遍会抑制潜艇尾流的流动涡结构，同时预报的宏观力亦存在明显的偏差；LES类湍流模型在捕捉潜艇流动涡结构、湍流应力方面表现优异，但是除LES-S模型外，针对宏观力的预测均出现明显的误差.因此在开展潜艇水动力及流场特征分析的综合研究中，LES-S模型是一个较为优异的选择.

4.3 潜艇周围涡场特征

基于上述分析得到的数值仿真策略对Sub-off的流场特征进行计算.由于艇体表面附体的存在，使得其周围流场较为复杂.图9所示为潜艇周围的涡系简图，涡系通过Q标准的重新采样体积进行表达，Q标准是Hunt等[19]提出的一项涡的判断方法.其中附着在围壳或尾舵附近的稍涡、马蹄涡(项链涡)会改变螺旋桨的来流均匀性而增加螺旋桨的载荷波动、提高空化风险及增加水动力噪声，湍流度巨大的涡结构会提高艇体表面激励力，引发艇身结构振动，削弱潜艇的隐身性能.艇体附着涡是由边界层流动引起的速度梯度变化形成，在边界层较为稳定的平行中体段强度微弱，在几何特征变化明显的首尾段有较强的艇体附着涡.图10所示为围壳附近多截面流场变量，图中：wi,wj,wk为涡量在3个方向的分量，图11对比围壳和尾舵表面的压力变化.从图10可以看出，围壳在轴向首先出现正压力梯度，在尾段出现负压力梯度，而其附近的涡量表明稍涡的脱落和压力梯度的变化存在密切关系，同时值得注意的是wk可能是稍涡变化的重要变量，存在明显脱落和分离的wi、wj是围壳随边涡变化产生的重要变量.由于尾舵和围壳相似的几何特征，本文认为尾舵稍涡和围壳稍涡的形成及影响因素相似，但尾舵稍涡的强度明显大于围壳稍涡.从图11可以看出，围壳和尾舵顶面和侧面的压力梯度差异可能是导致强度不同的关键因素， 需要注意的是图10和11中均对变量进行无因次化处理以方便对比分析.文献[17]的研究表明马蹄涡的存在与围壳(尾舵)、艇体边界层的三维流动有关，本文不对其进行深入研究.

图9 Sub-off周围涡系

图10 围壳附近多截面流场变量

图11 围壳和尾舵表面Cp分布

5 结语

本文以Sub-off为几何模型，针对特定工况，就10种湍流模型计算获得的结果进行分析，得出如下结论：

通过空间及时间收敛性分析，基于理查森外推法判定本文选用的网格方案具有良好的网格收敛性，最终采用了合适的网格开展后续的数值计算.本文选用的湍流模型中，LES-S模型在潜艇的宏观量如阻力、时均流场、伴流方面具有良好的性能或“随群性”；同时，在流场的二次量预报方面，LES类模型获得的结果较其他模型更加丰富.

潜艇涡系因其表面附体的存在较为复杂，但由于围壳和尾舵的几何相似性，其涡系具有内部形成机理的相似性，但因为附体装配位置及参数的差异，尾舵的涡系耦合和演化有较大差异.