0 引　言

操纵性是潜艇的重要性能之一，潜艇的操纵性对潜艇安全航行至关重要，因此研究潜艇操纵性一直是总体设计中的关键内容^[1]。潜艇斜航作为一种常见的影响潜艇操纵性的航行状态，受到了越来越多的学者关注。因此，对潜艇绕流场进行分析研究，探索其对水动力的影响，以此来提高潜艇操纵性有着极其重要的意义。

目前求取水动力系数的方法主要有3种^[2]：拘束或自由自航模型试验法、半理论半经验的估算方法、计算机数值模拟法。模型试验法是确定水动力系数最成熟的方法，测试结果被认为是精确度最高的数值通常用来作对比验证，但模型试验需要场地及操作人员同时耗费大量的时间和经济成本。半理论半经验的估算方法因采用经验公式大量节约了时间和经济成本，但其计算结果达不到工程上要求的精度。因此，在实际研究中常会将半理论半经验方法与其他方法结合运用，或者与其他计算方法进行精度对比，以确保得到更准确的结果。随着流体力学理论与计算机的进一步发展，需要大量计算的粘性流体数值模拟又被重视起来，并通过计算流体力学的CFD（Computational Fluid Dynamics）等算法，对潜艇操纵特性的数值模拟引起了更广泛的重视。郭海鹏等^[3]对SUBOFF潜艇在斜航工况下的操舵水动力及绕流场数值预报研究，研究了潜艇及尾舵不同舵角情况下受力与周围流场的变化。杜度等^[4]对均匀来流以及非均匀来流这2种来流方式下潜艇做直航运动以及斜航运动等4种工况下的周围流场进行数值模拟，探寻流场不均匀程度及漂角对潜艇受力的影响。

目前应用较多的RANS方法在小漂角情况下的预测结果已经能够达到一般的工程精度要求。但随着斜航漂角的增大使潜艇周围的流场变得异常复杂，发生流动分离现象，出现了涡的脱落，应用RANS方法目前无法捕捉分离涡系的演化过程^[5]。因此，本文采用DES和RANS两种不同的数值方法对潜艇绕流场进行分析研究，探索其对水动力的影响，为潜艇斜航大漂角情况下数值模拟提供参考。本文基于CFD软件STAR-CCM+采用RANS与DES两种数值方法对潜艇斜航工况下水动力及其流场进行研究。以Joubert BB2潜艇为研究对象，采用文献[6 - 7]里的试验工况对潜艇斜航进行数值模拟，探寻不同数值方法在不同舵角工况下对潜艇受力及周围流场变化的影响，为潜艇大漂角状态下的数值模拟以及数值方法的选取提供指导意见。

1 数值计算方法

本文采用基于雷诺平均（RANS）和分离涡模拟（DES）方法的求解器进行数值计算，并采用SST k-ω湍流模型实现方程的封闭。该湍流模型是一种综合了标准k-ω模型和k-ε模型优点的模型，具有以下特点^[8]。

1）保留了k-ω模型对近壁面处边界条件的良好处理能力：k-ω模型在边界层附近的流动处理方面表现良好，能够有效地捕捉壁面的湍流特性，提供准确的边界层预测。

2）避免了k-ω模型在自由剪切流区域对入口参数过于敏感的不足：相比于k-ω模型，SST k-ω模型能够更好地处理自由剪切流区域的湍流，降低对入口参数敏感度，提高计算的稳定性和准确性。

湍动能k和特定湍流耗散率ω的输运方程分别为：

$ \begin{gathered} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {U _i}k} \right)}}{{\partial {x_i}}} = {\widetilde Pk} - {\beta ^*}\rho k\omega + \frac{\partial }{{\partial {x_i}}} \left[ {\left( {\mu + {\sigma _k}{\mu _t}} \right)\frac{{\partial k}}{{\partial {x_i}}}} \right] ，\\ \\ \end{gathered} $

(1)

$ \begin{aligned}[b]\dfrac{{\partial}\left(\rho\omega\right)}{{\partial}t}+\dfrac{{\partial}\left(\rho U_i\omega\right)}{{\partial}x_i}= & \alpha\rho S^2-\beta\rho\omega^2+\dfrac{{\partial}}{{\partial}x_i}\left[\left(\mu+\sigma_{\omega}\mu_t\right)\dfrac{{\partial}\omega}{{\partial}x_i}\right]+ \\ & 2\left(1-F_1\right)\rho\sigma_{\omega2}\dfrac{1}{\omega}\dfrac{{\partial}k}{{\partial}x_i}\dfrac{{\partial}\omega}{{\partial}x_i}。\\[-12pt]\end{aligned} $

(2)

式中：$\rho $为流体密度；$ {U_i} $为速度矢量。其他未知量可参考相关文献[9]。

DES方法的SST模型只需要通过修改湍动能k的输运方程，将式(1)中的$ {\beta ^*}\rho k\omega $乘系数F_DES，得到新的输运方程，如下：

$ \begin{split} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {U_i}k} \right)}}{{\partial {x_i}}}=& {\widetilde P_k} - {\beta ^*}\rho k\omega {F_{\rm{DES}}} +\\& \frac{\partial }{{\partial {x_i}}}\left[ {\left( {\mu + {\sigma _k}{\mu _t}} \right)\frac{{\partial k}}{{\partial {x_i}}}} \right] ，\end{split} $

(3)

$ {{{F}}_{{\rm{DES}}}}{{ = {\mathrm{max}}}}\left[ {\frac{{{{{L}}_{{t}}}}}{{{{{C}}_{{\rm{DES}}}}{{\Delta }}}}{{(1 - }}{{{F}}_{{S}}}{{),1}}} \right]{{ }} 。$

(4)

式中：F_DES称为开关函数；C_DES为DES常数；$ \Delta= \sqrt[3]{\Delta x\Delta y\Delta z} $为网格长度尺度；F_S为混合函数，控制2种模型的切换。

2 数值计算设置 2.1 坐标系

本文采用2种右手直角坐标系对潜艇斜航进行数学描述，即固定在大地上的惯性坐标系O-xyz和坐标原点固定在艇体重心上的随艇一起运动的艇体坐标系O-x₁y₁z₁，如图1所示，其中u为来流速度。坐标系原点位于潜艇质心处，x轴正向由艇尾指向艇首，x₁轴正向指向潜艇运动方向；y轴、y₁轴正向指向潜艇右舷；z轴、z₁轴正向垂直向下。

2.2 研究对象

用于计算的潜艇模型为Joubert BB2潜艇（模型比例1∶52），该艇的模型参数如表1所示，三维几何形状如图2所示（由裸艇体，围壳以及4个尾舵X型对称布置在尾部）。

2.3 计算工况

具体计算工况参数如表2所示。考虑漂角（漂角的处理方法是将BB2以艇体局部坐标系O-xyz的Z轴为旋转轴进行旋转）在0°～–30°范围内的水动力性能，计算结果与试验^{[6 - 7]}数据进行对比，以此来验证数值方法的可靠性。

2.4 计算域划分

所选计算域及边界条件如图3所示。计算域总长为6L（L为潜艇总长，下同），左端边界距离艇首部1L设为速度进口，右端边界距离艇尾部4L设为压力出口，左右及上下边界距艇体均为1L均为速度进口。潜艇艇体表面均设定为无滑移壁面边界条件。

2.5 网格收敛性验证

在稳态计算过程中为消除网格质量和数量对计算结果带来的误差，需进行网格无关性验证^[10]。本文计算网格采用STAR-CCM+中的切割体网格、棱柱层网格以及表面重构的方法生成。在3个方向网格点布置依次根据$\sqrt 2 $倍变化的方式选取粗糙、中间及细密3种网格，3种网格的网格数量分别为896万、452万和237万。分别使用3种网格对潜艇直航状态下阻力进行数值计算。其计算结果如表3所示。

相邻的2套网格间阻力系数差用$ \varepsilon\mathit{_{\text{k}}} $表示，可以得到：

$ \varepsilon_{\mathit{\text{k}}21}=S_{\mathit{\text{k}}2}-S_{\mathit{\text{k}}1}=-0.016\times10^{-3}，$

(5)

$ \varepsilon_{\mathit{\text{k}}32}=S_{\mathit{\text{k}}3}-S_{\mathit{\text{k}}2}=-0.075\times10^{-3}\text{ }。$

(6)

由此可得收敛率R如下：

$ R=\dfrac{\varepsilon_{\mathit{\text{k}}21}}{\varepsilon_{\mathit{\text{k}}32}}=0.213。$

(7)

根据收敛准则0<R<1，单调收敛，满足网格收敛性要求。

3种网格计算结果误误差均在5%以内，网格数量增大使得计算结果更加精确，但考虑计算时间与计算精度二者之间的平衡，本文将采取S₂网格进行数值计算。

3 结果分析 3.1 水动力和力矩

为了便于观察将纵向力X、横向力Y和转首力矩N进行无因次化处理，无因次形式如下：

$ X' = X/\left(\dfrac{1}{2}\rho {V^2}{L^2}\right) ，$

(8)

$ Y' = Y/\left(\dfrac{1}{2}\rho {V^2}{L^2}\right)，$

(9)

$ N' = N/\left(\dfrac{1}{2}\rho {V^2}{L^3}\right)。$

(10)

式中：$\rho $为密度；V为速度；L为潜艇总长。

在图4中，可以看到不同漂角情况下潜艇所受纵向力X、横向力Y和转首力矩N的无因次化结果与试验值的对比情况。可以看出，漂角小于–15°时DES与RANS纵向力、横向力以及转首力矩系数预报精度均较好，误差均在10%以内。随着漂角的增大各个力和力矩系数的非线性增大，与试验结果误差逐渐增大尤其是受流动分离影响较大的纵向力。发生这种情况的原因之一是由于该试验使用的支撑塔以及整流罩对流场的不利影响^[7]，随着漂角的增大模型支撑结构负责在整流罩下游形成尾迹，这影响了模型后端的流场，该尾迹区域与艇体上的横流涡相互作用从而影响受力，使得模型试验的结果也存有相应的误差。通过2种方法数值计算结果的对比，总的来说DES方法的计算精度优于RANS方法。

3.2 压力分布

图5为RANS与DES方法不同漂角下潜艇表面压力分布情况，压力被表示为无因次的压力系数C_P，其无因次形式如下：

$ C_P=p/\left(\dfrac{1}{2}\rho V^2\right)。$

(11)

可以看出，斜航运动使得潜艇两侧的表面压力呈现出很强的非对称性，随着漂角的增大明显加剧了压力的非对称分布，即左右两侧的压力差增大，这就导致了潜艇的横向力以及转首力矩的增大。2种数值方法在相同舵角情况下的表面压力分布在主艇体部分大致相同，主要差异集中在潜艇围壳处。2种方法在漂角–20°时潜艇左侧围壳处压力分布差异最大，此时DES方法形成了一个低压区域阻挡了流体继续向下游流动从而使得潜艇在–20°时纵向力异常。总体来说，2种数值方法得到的表面压力分布差异较小，主要差异集中在潜艇围壳处。

3.3 流线分布

图6为RANS与DES方法在不同漂角情况下艇体周围（围壳及尾舵）的流线分布。可知，斜航运动绕流场呈明显的非对称性，不同漂角情况下的流场差异较大，同时同一漂角不同数值方法间流场也存在差异。RANS与DES方法在漂角为0°与–10°工况下，艇周围流线主要为附着流，2种方法流场几乎相同。在漂角为–15°工况下，RANS方法仍然是以附着流为主，DES方法捕捉到了轻微的流动分离，该情况并未对潜艇受力造成较大的影响。在漂角为–20°与–30°工况下，2种方法均发生了流动分离现象，其中DES方法的流动分离现象更加明显，这使得潜艇围壳及尾舵周围流速下降从而在其周围形成一定的涡，这些涡同时也会对潜艇产生相应的力，从而影响潜艇的受力，同时使得流向其下游区域的流体速度明显降低，阻碍了一部分本该流向尾舵部分的流体，从而影响尾舵处的流场。综上，随着漂角的增大潜艇周围流线分布越来越复杂，DES方法相较于RANS方法能够更好地捕捉到围壳及尾舵周围复杂的流场信息。

3.4 涡系结构

图7为不同漂角下船体周围的涡流结构，使用Q准则^[11]（Q=100 s⁻²）等值面进行展示。可以看出，在–10°漂角情况下主艇体的进流段和去流段均有一层厚厚的附着涡，指挥台围壳和尾翼与主艇体的接触部分产生马蹄涡。–20°及–30°漂角情况下围壳附近产生明显的涡脱落并向后延伸与尾舵处产生的涡相互作用从而形成复杂的涡系结构。漂角≤–10°情况下RANS与DES方法计算出的涡系结构差异不大，随着漂角的不断变大，RANS与DES方法得到的涡系结构差异越来越大，特是在大规模流动分离的情况下，RANS方法不能很好捕捉到围壳周围复杂的涡结构，而DES方法能够捕捉更为精细的涡结构。潜艇在粘性流体中运动时，流体会受到指挥台围壳和尾翼等障碍物的阻碍作用，导致流速降低，产生逆压梯度。当逆压梯度足够大时，流体的运动会被阻止，并且在逆压梯度的作用下，流体会向障碍物的上游方向运动，形成沿着围壳进流面根部周向分布的流动分离的涡旋^[12]。这种现象通常被描述为流动分离，是由于障碍物表面的逆压梯度导致流体分离并形成涡旋。这种流动分离会导致流体无法有效地进入围壳和尾翼的根部区域，从而影响潜艇的运动性能和流体动力学特性。综上，由于DES方法能够捕捉更为精细的涡结构所以其大漂角情况下的数值计算结果误差较小。

4 结　语

本文通过对Joubert BB2潜艇在斜航状态下不同漂角对潜艇周围的流场分布及其受力的影响进行分析，得出以下结论：

1）潜艇斜航工况下潜艇的力和力矩在水动力预报精度方面DES方法比RANS方法精度更高，DES方法的计算精度均在15%以内，RANS方法在大漂角情况下的计算精度超过了20%。

2）RANS方法和DES方法对潜艇斜航时围壳处压力分布的数值模拟结果差异相对不大，从数值结果来看DES方法优于RANS方法。

3）小漂角（漂角小于15°）情况下RANS与DES力和流场差距不大；大漂角（漂角大于15°）情况下RANS方法力和力矩的误差较大同时对流场的捕捉不够精细。DES方法对于流场捕捉更为精细，在水动力计算方面误差也更小，在潜艇复杂工况下的水动力机理研究具有明显优势。考虑到DES方法是建立在更高的时间成本上所计算出来的结果，因此建议在小漂角情况下可以采用RANS方法即可达到较好效果，随着漂角的增大RANS无法捕捉相应流场时采用DES方法进行计算这样能够在保证计算结果准确性的情况下节约一定的计算资源。