0 引　言

针对实尺度船舶的水动力性能研究，当前广泛采用的技术手段是利用模型尺度的水池试验或者数值模拟数据外推得到实船的水动力性能。然而，船模尺度与实船尺度之间往往不满足某些相似准则，其中最主要的是表征流体惯性力与粘性力比值的雷诺数。对于不同尺度下的船模试验或数值计算，不能同时满足傅汝德数与雷诺数相等，因此在水动力外推预报时存在“尺度效应”^[1]。对于水下航行体而言，“尺度效应”问题是影响其水动力性能预报精度的重要因素，给实尺度船舶航行性能评估带来一定困扰。随着计算流体动力学（computational fluid dynamics, CFD）的发展，应用CFD方法对船舶及水下航行体水动力的尺度效应分析并预报实船水动力性能是当前船舶水动力学研究领域的热点之一。

操盛文等^[2]对高雷诺数下SUBOFF潜艇的全附体模型绕流进行数值模拟，研究了网格数量及艇体尺度对高雷诺数条件下潜艇阻力计算结果的影响。吴方良等^[3]对SUBOFF潜艇的主艇体模型和全附体模型在不同雷诺数条件下的流场进行了数值模拟，分析了潜艇粘压阻力系数随雷诺数变化的规律。司朝善等^[4]以SUBOFF潜艇全附体模型为对象，通过改变潜艇尺度实现不同雷诺数下的潜艇直航绕流数值模拟，分析了模型尺度对潜艇阻力性能和绕流场的影响规律。王展智等^[5]以DTMB5415裸船体和全附体模型为研究对象，计算了包含实尺度在内的多种缩尺模型的流场，详细分析了双桨水面舰船附体阻力的尺度效应。张恒等^[6]以KCS船型作为研究对象，计算了不同雷诺数下的阻力并对船体尾部流场进行分析，探讨了改变介质的粘性系数以满足全相似在实尺度预报方面的可行性。蔡博奥等^[7]针对不同缩尺比的三体船模型进行了各阻力成分的数值预报，探讨了三体船阻力成分关于水动力系数的变化规律并提出了三体船阻力换算的改进方法。师超等^[8]开展了不同船体尺度下的裸船体斜航和定常回转数值计算，探讨了不同尺度下船体操纵运动水动力的变化规律，并结合船舶操纵运动数学模型分析了尺度效应对操纵性的影响。Dogrul等^[9]对不同尺度下的KCS和KVLCC2船的阻力试验进行了数值模拟，探讨了不同阻力成分的尺度效应问题，并据此分析了不同外推方法在船舶实尺度阻力预报方面的有效性。宋科委等^[10]以DTMB5415船型为研究对象，预报了多种尺度、多种工况及多种船体壁面条件下的阻力性能，对船体形状因子、兴波阻力及摩擦阻力的尺度效应问题进行了分析，探讨了多种实船阻力性能预报方法存在差异的原因。Sezen等^[11]以SUBOFF潜艇模型为对象，研究了尺度效应对潜艇阻力等水动力性能参数的影响，通过对实尺度潜艇阻力试验的数值模拟，验证了基于傅汝德外推法的实尺度阻力预报方面的有效性。

从以上研究可以看出，针对船舶及水下航行体水动力性能尺度效应的研究大多集中在阻力性能方面，而针对船舶及水下航行体操纵性能尺度效应问题的研究还十分有限。本文针对水下航行体的操纵性能评估需求，开展水下航行体操舵水动力性能的尺度效应问题数值研究，应用CFD方法开展不同尺度的水下航行体直航操舵试验数值模拟，探讨艇体绕流场及水动力特性随雷诺数的变化规律，为构建实尺度水下航行体操舵水动力性能的外推预报方法提供参考。

1 数值方法

考虑到水下航行体操舵工况下存在的大尺度流动分离现象，应用分离涡模拟（Detached Eddy Simulation, DES）求解器开展数值计算。求解器所采用的DES方法基于剪切应力输运（Shear Stress Transfer, SST） $ k - \omega $ 湍流模型进行构造。SST $ k - \omega $ 两方程湍流模型如下：

$ \left\{ \begin{gathered} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_j}k} \right)}}{{\partial {x_j}}} = P - {\beta ^*}\rho \omega k + \\ \qquad\qquad\qquad\qquad\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + {\sigma _k}{\mu _t}} \right)\frac{{\partial k}}{{\partial {x_j}}}}\right]，\\ \frac{{\partial \left( {\rho \omega } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_j}\omega } \right)}}{{\partial {x_j}}} = \frac{{\rho \gamma }}{{{\mu _t}}}P - {\beta ^*}\rho {\omega ^2} + \\ \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + {\sigma _\omega }{\mu _t}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + 2\left( {1 - {F_1}} \right)\frac{{\rho {\sigma _{\omega 2}}}}{\omega } \cdot \frac{{\partial k}}{{\partial {x_j}}} \cdot \frac{{\partial \omega }}{{\partial {x_j}}} 。\\ \end{gathered} \right. $

(1)

式中： $\ \rho $ 为流体密度； $ {u_j} $ 为速度矢量； $ {x_j} $ 为位置矢量； $ k $ 和 $ \omega $ 分别为湍流动能和湍流比耗散率； $\ \mu $ 和 $ {\ \mu _t} $ 分别为层流和湍流黏性系数。源项 $ P $ 和函数 $ {F_1} $ 的具体形式以及系数 $ {\ \beta ^*} $ ， $ \gamma $ ， $ {\sigma _k} $ ， $ {\sigma _\omega } $ 和 $ {\sigma _{\omega 2}} $ 参考文献[12]。

将RANS湍流模型中的长度尺度采用DES的长度尺度代替，即可获得基于该湍流模型的DES方法，DES的长度尺度表示如下：

$ {l_{DES}} = \min \left( {{l_{RANS}},{l_{LES}}} \right)，$

(2)

其中，RANS长度尺度为最近壁面距离，表达如下：

$ {l_{RANS}} = \frac{{{k^{1/2}}}}{{{\beta ^*}\omega }}，$

(3)

此外，LES长度尺度仅与网格间距有关，表达如下：

$ {l_{LES}} = {C_{DES}}\Delta 。$

(4)

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

由DES模型的构造可知，当RANS长度尺度激活时，以RANS湍流模型的方式计算；而当LES长度尺度激活时，湍流模型转换为经典的Smagorisnky亚格子模型，从而在求解流场时实现RANS方法与LES方法的耦合。

2 研究对象及工况

本文以美国DARPA潜艇模型SUBOFF-8为研究对象，潜艇几何形状和主要尺度如图1和表1所示。针对不同缩尺比的潜艇模型，开展直航工况下的操舵试验数值模拟，工况见表1。数值模拟工况满足傅汝德数相似，对应傅汝德数Fr=0.512。

主要关注对水下航行体操纵性能有显著影响的纵向力 $ X $ 、横向力 $ Y $ 及转首力矩 $ N $ 。为便于对比分析，水动力及力矩均采用无因次形式描述，表达如下：

$ {X^{'}} = \frac{X}{{0.5\rho {V^2}{L^2}}} ，$

(5)

$ {Y^{'}} = \frac{Y}{{0.5\rho {V^2}{L^2}}}，$

(6)

$ {N^{'}} = \frac{N}{{0.5\rho {V^2}{L^3}}} 。$

(7)

3 计算域与网格生成

计算域为长方体，如图2所示。坐标系原点位于潜艇尾端，计算域入口（Inlet）边界距坐标系原点2.5倍艇长L_OA，计算域出口（Outlet）边界距坐标系原点4.5倍艇长L_OA，计算域左右边界（Left、Right）与上下边界（Top、Bottom）距坐标系原点2倍艇长L_OA。边界条件方面，出口边界设为压力出口，其他边界均设置为速度入口，艇体表面设为无滑移壁面。

计算域采用非结构六面体网格、棱柱层网格以及表面重构的方法进行网格划分。为了保证潜艇附近的网格质量，对其表面进行了加密。通过使用切割体网格单元，在艇体附近及其尾流进行网格密化，从而更好地捕捉流动细节，获得更准确的流场信息。通过调整棱柱网格层的厚度、层数及增长率来控制近壁面第一层网格高度，使无因次壁面距离y+值在30～300之间。图3和图4给出了潜艇表面和计算域的网格分布情况。

4 网格收敛性分析

采用Stern等^[13]提出的方法开展网格收敛性分析：

$ {\varepsilon _{G21}} = {S_{G2}} - {S_{G1}}，$

(8)

$ {\varepsilon _{G32}} = {S_{G3}} - {S_{G2}}，$

(9)

$ {R_G} = \frac{{{\varepsilon _{G21}}}}{{{\varepsilon _{G32}}}}。$

(10)

式中： $ {S_G} $ 为该网格下的参数计算结果； $ {\varepsilon _G} $ 为相邻2套网格对应参数值的差； $ {R_G} $ 为网格收敛率。

收敛性分析结果存在3种情况：1）0< $ {R_G} $ <1，单调收敛；2） $ {R_G} $ <0，振荡收敛；3） $ {R_G} $ >1，网格发散。当网格发散时，无法估计不确定性。网格振荡收敛时，不确定性由下式估计：

$ {U_G} = \left| {\frac{1}{2}\left( {{S_U} - S{}_L} \right)} \right|。$

(11)

式中： $ {S_U} $ ， $ S{}_L $ 分别为最密网格计算过程中的上、下限； $ {U_G} $ 为网格不确定性。

网格收敛时，采用广义理查森外推法计算误差的单项估算值和准确度阶数：

$ \delta _{R{E_{G1}}}^{*\left( 1 \right)} = \frac{{{\varepsilon _{{G_{21}}}}}}{{r_G^{{p_G}} - 1}} ，$

(12)

$ {p_G} = \dfrac{{\ln \left( {\dfrac{{{\varepsilon _{{G_{32}}}}}}{{{\varepsilon _{{G_{21}}}}}}} \right)}}{{\ln \left( {{r_G}} \right)}} 。$

(13)

式中， $ {r_G} $ 为网格细化比。

若误差的单项估算值 $ \delta _{R{E_{G1}}}^{*\left( 1 \right)} $ 和准确度阶数 $ {p_G} $ 已知，可以估计不确定性，修正因子定义为：

$ {C_G} = \frac{{r_G^{{p_G}} - 1}}{{r_G^{{p_{Gest}}} - 1}}，$

(14)

式中， $ {p_{Gest}} $ 为所采用数值方法的极限或理论精度，一般可以选择离散精度阶数。

如果修正因子远离1，缺少置信度，数值不确定度可以定义为：

$ {U_G} = \left\{ {\begin{array}{*{20}{c}} {\left( {9.6{{\left( {1 - {C_G}} \right)}^2} + 1.1} \right)\left| {\delta _{R{E_{G1}}}^*} \right|,\left| {1 - {C_G}} \right| < 0.125}，\\ {\left( {2\left| {1 - {C_G}} \right| + 1} \right)\left| {\delta _{R{E_{G1}}}^*} \right|,\left| {1 - {C_G}} \right| \geqslant 0.125} 。\end{array}} \right. $

(15)

在进行网格收敛性分析时，采用粗（S₃）、中（S₂）、细（S₁）三套网格，网格细化比 $ {r_G} $ = $ \sqrt 2 $ ，不同缩尺比下的网格数量如表2所示。不同缩尺比下 $ X $ ， $ Y $ 和 $ N $ 的网格不确定性分析如表3～表7所示。由表可知，当L/L_OA=1时， $ X $ ， $ Y $ 和 $ N $ 的网格收敛率均满足 $ {R_G} $ <0，即网格振荡收敛，不确定度为2.83%～8.48%；当L/L_OA=2时， $ X $ ， $ Y $ 和 $ N $ 的网格收敛率均满足0< $ {R_G} $ <1，即网格单调收敛，准确度阶数为2.40～4.51，不确定度为1.54%～5.53%；当L/L_OA=4时， $ X $ ， $ Y $ 和 $ N $ 的网格收敛率均满足 $ {R_G} $ <0，即网格振荡收敛，不确定度为3.12%～4.54%；当L/L_OA=8时， $ X $ 的网格收敛率 $ {R_G} $ <0，即网格振荡收敛，不确定度为4.26%，但是 $ Y $ 和 $ N $ 的网格收敛率 $ {R_G} $ >1，即网格发散；当L/L_OA=16时， $ X $ 的网格收敛率 $ {R_G} $ <0，即网格振荡收敛，不确定度为4.09%，但是 $ Y $ 和 $ N $ 的网格收敛率 $ {R_G} $ >1，即网格发散。

对于网格收敛性研究中存在网格发散的情况，其原因是多方面的：潜艇大舵角操舵状态下舵附近流场往往存在大尺度流动分离现象，流场具有显著的非线性和非定常性；所采用的DES方法是一种混合算法，其湍流建模与网格尺度存在关联性，不同密度网格其流场内局部流动的湍流求解方式也有所差异，导致难以呈现出较好的网格收敛性。

5 计算结果 5.1 泄出涡系结构

图5～图9给出了不同缩尺比潜艇的艇体及舵附近的无量纲Q准则瞬时等值面。等值面取Q=100s⁻²，并采用无因次纵向速度U_x/U₀着色。可以看出，垂直舵附近产生了明显的泄出涡结构，这主要是由于舵背风面处由于舵角所产生的流动分离现象所致。通过观察不同缩尺比的等值面图可以发现，舵附近的泄出涡规模随着缩尺比的改变而逐渐变化。在缩尺比较小（L/L_OA=1）时，舵背风面产生了大规模的流动分离现象，泄出涡脱落起始位置接近舵导缘位置。同时，舵叶梢部泄出涡与背风面泄出涡发生明显的相互作用，导致梢部泄出涡结构不明显。随着缩尺比增大（L/L_OA=2，4，8，16），泄出涡脱落起始位置向尾缘移动，背风面泄出涡规模逐渐减小。同时，舵叶梢部泄出涡结构逐渐清晰，受背风面泄出涡的影响逐渐减小。总体而言，在舵角不变的条件下，舵背风面流动分离现象随着雷诺数的增大逐渐减弱。

5.2 舵面压力分布

图10～图14给出了不同缩尺比潜艇的舵中剖面处的压力分布。压力采用无因次系数 $ {C_P} $ 表示， $ {C_P} = P/ $ $\left( {0.5\rho {V^2}} \right) $ 。横坐标用无因次量x/C表示，x/C=0表示舵前缘位置，x/C=1表示舵尾缘位置。在靠近舵前缘位置，舵面分别达到正压和负压峰值。在舵剖面中后区域压力分布趋于平缓。从图中可以看出，不同缩尺比下舵剖面压力分布存在显著差异，特别是舵背风面。在L/L_OA=1时，舵背风面压力分布以x/C=0.2为分界线呈现出不同的分布，x/C=0～0.2区域，表现为负压峰值；x/C=0.2～1.0区域，舵面压力趋近于0，且分布十分平缓。随着L/L_OA逐渐增加，舵背风面压力分布开始发生变化，分界线逐步向尾缘移动，且靠近尾缘部分舵面压力开始出现不规则波动，但整体压力差有所增大。

5.3 水动力及力矩

图15～图17给出了作用在不同缩尺比潜艇上的无因次纵向力、横向力及转首力矩的计算结果。对于纵向力，在雷诺数Re=1.4×10⁷～31.3×10⁷范围内变化最为明显，这意味着该雷诺数区间内尺度效应的影响十分显著。随着雷诺数进一步增加，纵向力系数的变化趋缓，尺度效应的影响有所下降，但仍然难以忽略。对于横向力，在雷诺数Re=1.4×10⁷～3.9×10⁷范围内变化较为明显，随着雷诺数进一步增加，横向力系数的变化显著减缓，意味着尺度效应的影响明显减弱。此外，转首力矩与横向力的变化规律基本一致。从以上分析可以看出，与潜艇阻力性能相关的纵向力系数其尺度效应较为明显，且随雷诺数的增加其影响变化较小。相较而言，与潜艇操纵性能相关的横向力及转首力矩系数其尺度效应在雷诺数高于某一值后影响显著减弱，横向力及转首力矩系数趋近于某一定值，这也意味着在研究潜艇操纵运动特性时，尺度效应的影响不可忽略。

6 结　语

本文以SUBOFF潜艇模型为研究对象，按照1，2，4，8，16的缩尺比对潜艇进行等比例缩放，采用CFD方法对不同尺度的潜艇模型的直航操舵试验进行数值模拟，对获得的绕流场和水动力特性进行分析，可得到以下结论：

1）不同雷诺数下的艇体附近流场存在显著差异。舵面附近的流动分离现象随雷诺数的增大而减弱，流动分离发生位置逐步后移，导致泄出涡系结构也有所减弱，其中雷诺数较小时变化更为明显。舵面压力分布也随之发生改变，舵背风面的负压区逐渐扩大，导致舵剖面压差随雷诺数增大而增大。

2）潜艇纵向力、横向力及转首力矩系数随雷诺数的变化而发生显著变化。随雷诺数增大，尺度效应对潜艇纵向力的影响虽有所下降但依然显著，而对横向力和转首力矩的影响显著减弱。