潜艇斜航工况下操舵水动力及绕流场数值研究

郭海鹏，祁江涛，袁文鑫，李广年

(1.宁波大学 海运学院，浙江 宁波 315211；2.中国船舶科学研究中心，江苏 无锡 214082)

0 引 言

操纵性是衡量潜艇综合性能优良的重要指标之一，良好的操纵性能是潜艇安全航行的重要保证[1-2]。潜艇操纵工况下绕流特性极为复杂，不但有主附体结合部和流水孔处各种形式涡的产生与演化，而且还涉及艇体、舵以及螺旋桨之间的复杂相互干扰，导致针对潜艇操纵性的研究具有很高的难度。早期研究主要依靠模型试验方式进行，对人力物力要求高而且试验设施复杂，难以在潜艇设计阶段进行广泛应用。

随着计算机技术的发展，计算流体动力学（computational fluid dynamics, CFD）方法在水下航行体操纵性预报中发挥越来越重要的作用。张楠等[3-5]提出了适用于水下航行器流场预报的数值方法，并开展了带自由液面、全附体模型的数值模拟。孙铭泽等[6]实现了全附体潜艇平面运动机构试验数值模拟，并将计算得到的水动力导数与试验值以及直线拖曳试验或回转试验的数值计算结果进行了对比验证。庞永杰[7]以椭球体为研究对象计算了其水动力导数，并通过与理论值进行对比验证了所采用数值方法的有效性。焦玉超等[8]对不同X 形尾舵的潜艇绕流场进行了数值模拟，分析了X 舵角度对潜艇性能的影响。翟朔等[9]对某潜艇模型的尾部水平操纵面分别进行了共翼型设计和非共翼型设计，并计算了2 种操纵面产生的艇体水动力和尾流特征。

从目前已有研究可以看出，应用CFD 方法开展潜艇操纵性能预报及评估已经成为领域内的趋势和研究热点。本文基于CFD 平台STAR-CCM+对潜艇斜航工况下的操舵水动力及绕流场进行数值研究，以国际上广泛用于对比验证研究的SUBOFF 潜艇模型为研究对象，分别应用雷诺平均（reynolds-averaged navier-stokes,RANS）与分离涡模拟（detached eddy simulation, DES）2 种方法对潜艇不同舵角下的斜拖试验进行数值模拟，据此分析潜艇纵向力、横向力、转首力矩以及绕流场随舵角的变化规律，揭示潜艇操纵工况下水动力特性与绕流场间的内在联系，为水下航行体大舵角状态操纵性能预报和评估提供技术支持。

1 数值计算方法

采用STAR-CCM+中的RANS 与DES 求解器开展相关数值计算工作，其中RANS 方法和DES 方法均采用剪切应力传输模型（shear stress transfer, SST）k-ω湍流模型进行构造。SSTk-ω两方程湍流模型为：

式中：ρ为流体密度；uj为速度矢量；xj为位置矢量；k和 ω分别为湍流动能和湍流比耗散率；μ和 μt分别为层流和湍流黏性系数。源项P和函数F1的具体形式以及系数 β*，γ，σk，σω和σω2可参考SSTk-ω模型的相关文献[10]。

将RANS 湍流模型中的长度尺度使用DES 的长度尺度代替，就获得了基于该湍流模型的混合RANS/LES方法即DES 方法。DES 的长度尺度由RANS 和LES 的长度尺度以如下方式混合得到，即

LES 长度尺度由网格间距得到，即

式中，Δ为局部网格间距，一般取网格单元与相邻网格单元中心连线长度的最大值。

RANS 长度尺度为最近壁面距离，即

由DES 模型的构造可知，当RANS 长度尺度激活时，湍流模型为传统的RANS 湍流模型，而当LES 长度尺度激活时，湍流模型转换为经典的Smagorisnky 亚格子模型，从而实现RANS 和LES 方法的转换。

2 数值计算设置

2.1 研究对象

以美国DARPA 潜艇模型SUBOFF-8 为研究对象，该潜艇由裸艇、指挥台围壳及4 个完全相同的尾翼（左右为升降舵，上下为方向舵）组成，主艇体总长Loa=4.356 m，垂线间长L=4.261 m，其中前体长1.016 m，平行中体长2.229 m，后体长1.111 m；艇体最大直径D=0.508 m；指挥台围壳长0.368 m，高0.205 m，上部有2∶1 的椭圆形横截面的顶盖；4 个相同的尾翼截面为NACA0020，呈十字布置在尾部，顶端弦长为0.152 m，详细几何特征由文献[11]提供，三维几何形状如图1所示。

图1 SUBOFF 几何模型Fig.1 Geometric model of SUBOFF

2.2 计算工况

计算工况具体参数如表1 所示。考虑潜艇在漂角β=-4°、舵角0°～40°范围内的水动力性能，其中δ在0°～15°范围内数值计算结果可与Roddy[12]所做的公开实验数据进行对比，从而确定数值方法的可靠性。在STAR-CCM+中，漂角的处理方法是将SUBOFF 以艇体质心所在的Z轴为旋转轴，进行逆时针旋转4°，其他保持不变。

表1 工况设定Tab.1 Setup of study condition

2.3 计算域设置

计算域总长为7 倍艇长，首部距离进口1.5 倍，尾部距离出口4.5 倍，左右及上下边界距艇体均为2.0 倍艇长。对于边界条件的设定，除计算域出流面采用压力出口条件以外，其他边界均采用速度入口条件，艇体所有表面采用无滑移壁面边界条件。计算域尺度及边界条件如图2 所示。

图2 计算域尺度及边界条件Fig.2 Dimensions and boundary conditions of the computational domain

计算域采用STAR-CCM+中的切割体网格进行离散。切割体网格属于非结构网格，适用于具有复杂外形模型的网格划分。切割体网格在生成时可以通过控制体进行任意区域的网格加密，同时还可以对几何模型中的控制线、控制面进行线加密及面加密。为了保证潜艇附近的网格质量，对其表面进行了加密；通过使用切割体网格单元，在艇体附近及其尾流进行网格密化，从而更好地捕捉流动的细节，获得更准确的流场信息。同时，在近壁面生成棱柱层网格以便于边界层流动的求解。通过调整棱柱层网格的总厚度、层数及增长比率来控制近壁面第一层网格高度，使近壁面Y+值在30～60 范围内。网格数量在500 万左右，网格划分情况如图3 所示。

图3 计算域网格划分情况Fig.3 Gird generation in the computational domain

3 结果分析与讨论

3.1 水动力及力矩

图4 给出了RANS 与DES 方法获得的作用在潜艇上的纵向力X、横向力Y及转首力矩N。为了便于分析，均采用无因次形式表述：

图4 数值结果与模型试验数据对比Fig.4 Comparison between the numerical results and the experimental data

式中：ρ为流体密度；V为航行速度；L为潜艇总长。

相对X′而言，无论RANS 还是DES 对阻力的预报精度都在±3%以内；δ的变化会对阻力产生较大的影响，δ在0°～30°范围内，随着舵角的增大，艇体所受的阻力以抛物线的趋势增大；随着δ的进一步增大，RANS 计算结果在30°～40°范围内基本呈线性增大，DES 计算结果先有略微减小然后再增大。相对Y′而言，RANS 计算结果与试验值的吻合程度略好于DES 方法，RANS 计算的误差均在10%以内，DES 计算除了δ=15°时误差偏大，其余工况均在12% 以内；δ在0°～30°范围内，DES 与RANS 计算得到的横向力随着舵角的增大而近似呈现线性增大；δ在30°～35°范围内，RANS 及DES 计算的横向力显著下降，发生了明显的舵失速现象；进一步增大舵角，艇体所受的横向力有微幅的增加，这表明舵效并未随舵角进一步增大而发生明显改善。相对于N′而言，DES 计算在各工况下的计算精度均在±12% 以内，而RANS 方法在δ=10°和δ=15°计算误差较大，均超过20%；N′与Y′的变化规律保持基本一致。通过对水动力分析可知，就对比现有的公开试验数据而言，RANS 方法在X′与Y′相关的计算精度比DES 方法略好，但是DES 方法对N′的计算精度明显优于RANS 方法。

3.2 泄出涡系结构

图5 给出了由RANS 及DES 方法获得的不同舵角下艇尾区域的涡系结构，其中涡系结构采用Q准则等值面表示，并以无量纲纵向速度Ux/U0着色。可以看出，随着方向舵舵角的变化，艇尾区域涡系发展呈现出显著的差异。在小舵角（δ=10°）下，舵背风面主要以附着流为主，舵梢部有轻微涡泄出。随着舵角增大（δ=20°），方向舵附近产生的涡系结构明显增多。方向舵背风面靠近尾部区域开始产生泄出涡系。同时，由于偏转方向舵的阻滞效应，在艇身与方向舵的交汇处形成了明显的马蹄涡结构。随着舵角的进一步增大（δ=30°），舵背风面完全由泄出涡系覆盖，而且马蹄涡结构也更为显著。当舵角达到δ=40°时，舵背风面的泄出涡完全脱离舵表面发展，在下游形成了复杂的涡系结构，与舵失速状态相对应。通过对比可知，在小舵角工况下，RANS 方法与DES 方法得到的艇尾区域涡系结构差异不明显，随着舵角增大，二者的差异越发明显。特别是发生大规模流动分离的工况下，DES方法捕捉到的涡系结构更加精细。

图5 不同舵角工况下的Q=100 s-2 等值面图Fig.5 Iso-surface of Q=100 s-2 under different rudder angles

3.3 流线分布

图6 给出了DES 与RANS 方法获得的舵叶中间位置处的舵附近流线分布。可知，在舵角δ=10°工况下，舵剖面附近流动主要以附着流为主，而且2 种方法得到的流场没有明显的差异。在舵角δ=20°工况下，舵剖面背流面发生了明显的流动分离现象，其中RANS方法呈现出的流动分离更为明显，导致该区域的速度明显下降。在舵角δ=30°工况下，舵剖面的整个背流面均发生了流动分离，2 种数值方法获得的流线分布差异也更加明显。RANS 方法中背流面处的流动分离更加明显，而DES 方法捕捉到的流动分离并未随舵角明显增强。另外，2 种方法捕捉得到流场均在舵下游呈现出明显的低速区。在舵角δ=40°工况下，舵剖面背流面处的流动分离显著增强，对附近流场产生的扰动也十分明显，在下游形成了更为显著的低速区，这一流场特征与舵失速工况下的阻力增加及升力下降相对应。综上，DES 方法在流场捕捉方面具有明显优势，能清晰地捕捉艇尾区域的复杂流场信息。

图6 不同舵角工况下舵附近流线分布Fig.6 Streamlines around the rudder under different rudder angles

3.4 压力分布

图7 给出了RANS 方法与DES 方法获得的不同舵角工况下的舵剖面压力分布情况。压力采用无因次系数Cp表示，Cp=p/(0.5ρV2)。横坐标由无因次量x/C表示，坐标原点位于该舵前缘处，无因次量“1”表示舵尾缘。在舵角δ=10°工况下，舵剖面压力呈现出明显的吸力面和压力面，正压和负压的峰值均位于前缘附近，而且2 种方法获得的压力分布较为一致。在舵角δ=20°工况下，DES 方法得到的吸力面负压区较RANS方法更为显著，但整体分布差异很小。在舵角δ=30°工况下，吸力面负压区更为显著，同时发生流动分离的区域压力分布较为平缓。DES 方法得到的吸力面负压区明显强于RANS 方法，这也与该舵角工况下DES计算得到的舵横向力高于RANS 计算结果相对应。在舵角δ=40°工况下，吸力面负压区峰值明显下降，而且发生流动分离的区域压力分布呈平直状，导致吸力面与压力面之间的压力差显著下降，这与舵发生失速导致的舵效下降相对应。总体而言，2 种数值方法得到的压力分布差异相对较小，与潜艇操舵水动力性能的一般规律较为吻合。

图7 不同舵角工况下舵剖面压力分布图Fig.7 Pressure distribution on the rudder section under different rudder angles

4 结 语

本文针对潜艇斜航工况下的操舵水动力性能及绕流场进行数值模拟研究，主要得到以下结论：

1）潜艇斜航工况下的纵向力、横向力及转首力矩随着舵角的变化呈现复杂的变化规律，在潜艇操舵水动力预报精度方面DES 方法整体上优于RANS 方法；

2）DES 方法对于潜艇操舵工况下艇尾区域的绕流场捕捉更为精细，在潜艇操纵工况下的水动力学机理探讨方面有明显的优势，在潜艇复杂操纵运动分析和性能评估方面具有很好的工程应用潜力。