黄永生，杨晨俊，董小倩

1 上海交通大学 海洋工程国家重点实验室，上海 200240

2 高新船舶与深海开发装备协同创新中心，上海 200240

0 引 言

快速性预报是水面船舶及水下航行体设计的重要环节之一，其传统做法是在阻力、敞水及自航模型试验的基础上对阻力、伴流及敞水性能的尺度效应采用1978 ITTC规程等经验方法进行修正。对于常规螺旋桨推进的船舶，采用该方法虽然可靠性好，但成本较高。高速水下航行体一般采用对转桨（CRP）或泵喷推进，推进器部件间的相互作用给尺度效应修正带来了一定的困难（主要是伴流）。另一方面，随着计算流体力学（CFD）的快速发展，粘流CFD计算技术在快速性预报方面的应用日益增多，并有逐步取代模型试验的趋势。目前，阻力和敞水性能的数值计算精度已基本达到工程实用要求，自航计算的精度也在逐步提高。与模型试验相比，粘流CFD计算技术的另一个优势在于能够进行实尺度模拟，从而无需进行尺度效应修正。实尺度模拟对计算机容量和速度的要求较高，其主要困难在于船体，而推进器则相对容易实现。

自航计算需要同时模拟船体与推进器，网格量和计算量均较大；从理论上看，需要进行非定常模拟，但所需计算时间会更长。因此，一般采用定常或准定常的方法模拟船体与推进器的相互作用。为了降低网格生成难度，节省计算时间，进行定常计算时可以采用体积力模型来代替螺旋桨的作用。体积力模型分描述型和迭代型2种，前者的体积力径向分布形式是基于经验公式，后者则是通过螺旋桨升力面、面元等势流方法计算得到。基于描述型体积力模型，Choi等［1］对VLCC进行了自航模拟和快速性预报；傅慧萍［2］针对KCS进行了带自由液面的自航模拟，研究了扭矩对计算精度的影响；吕晓军等［3］开展了潜艇自航计算研究。Kim等［4］的研究表明，由于迭代型体积力模型是根据桨叶的真实几何及实效伴流分布来计算体积力的分布，所以它对艉部流场的模拟比描述型体积力模型的结果更加准确。虽然采用体积力方法对自航因子的预报结果较好，但毕竟是用力场代替了真实桨叶，对于桨叶对周围流体的排挤效应和扰动难以精确体现，从而限制了计算精度和流场细节。2005年以后，采用螺旋桨真实几何的数值自航模拟方法逐步发展起来。Lübke［5］应用CFX软件对KCS进行了自航模拟；Choi等［6］对一系列船型进行数值自航模拟和快速性预报，充分验证了其计算方法的精度和可靠性；熊鹰等［7］则针对傅汝德数较低的排水型船舶提出了一种简化的自由面处理方法，用以加快计算速度；杨琴等［8］针对配有七叶桨的SUBOFF潜艇进行了自航模拟，并分析给出了自航因子。上述自航模拟均是在模型尺度下进行，仍需通过修正尺度效应来预报实船快速性。随着计算机硬件能力的提高，近年来，有关实尺度自航计算的研究逐步增多。Castro等［9］以KCS为对象开展了实尺度自航模拟，并研究了推进因子的尺度效应；李亮等［10］也对KCS开展了类似的研究。

目前，针对高速水下航行体推进性能方面的数值计算研究相对较少。何文生［11］较早地采用螺旋桨的真实几何对高速水下航行体进行了RANS模拟，但受计算机能力的限制，简化了附体。张涛等［12-13］针对对转桨的敞水性能，开展了RANS计算研究。但是，针对高速水下航行体的自航计算研究，在模型和实尺度方面均未见公开发表的文献。

为此，本文拟开展对转桨推进的高速水下航行体实尺度自航CFD计算研究，基于阻力、敞水及自航模拟结果，分析自航因子。将计算处于深潜状态的航行体，而不考虑自由液面和空化的影响。

1 计算方法

本文采用求解RANS方程的方法，来开展对转桨推进的高速水下航行体自航性能实尺度计算研究。假定流体不可压缩，流动为全湍流，采用标准k-ε湍流模型和非平衡壁面函数。控制方程包括连续性方程、动量输运方程及湍流输运方程，此处不再赘述。网格生成采用ICEMCFD及GAM⁃BIT软件，数值计算采用FLUENT软件。

1.1 计算域及其离散

如图1所示（图中D0为航行体平行中体段的直径；L为航行体长度），阻力及自航计算域为包围航行体的回转体，但阻力计算域不包含前、后桨叶。将计算域划分为4个子域：Fore-body和Rudder/fin为静止域；Rotor-F和Rotor-A在阻力计算中为静止域，在自航计算中为旋向相反的旋转域。4个子域之间依次设有Int.mid，Int.aft和Int.CRP三对交界面。

子域Fore-body和Rudder/fin采用结构化单元离散，对航行体头部及附体等几何变化剧烈的区域进行了网格加密。用圆柱面CS将子域Rotor-F和Rotor-A分割为2个部分，其中CS的外部用结构化单元离散；用通过尾端面的横截面XS进一步分割子域Rotor-A，其中CS内部和XS下游部分用棱柱单元离散。在自航计算中，CS内部以及包含前、后桨叶的区域采用四面体及五面体非结构化单元离散；在阻力计算中，原前、后桨叶所在区域用结构化单元离散。图2所示为阻力、自航计算关键区域的物面/截面部分区域面网格。

如图3所示，对转桨的敞水计算采用圆柱形全流道计算域，图中，DF为前桨直径。前、后桨之间的交界面Int.CRP将计算域分为了2个旋转方向相反的子域，每个旋转域又进一步划分为3个子域。其中子域Far-F和Far-A采用六面体结构化单元离散；包围前桨及毂帽的子域Rotor-F和包围后桨的子域Rotor-A采用四面体及五面体非结构化单元离散；子域UpS和DownS采用棱柱体单元离散。为了保证桨叶间网格的周期性，对于所有子域，在相邻桨叶之间生成周期面，将全流道沿周向等分成Z（Z为桨叶数）个子流道，通过复制子流道网格、合并周期面网格，得到全流道网格。图4所示为子域UpS，Rotor-F，Rotor-A及DownS的子流道几何。

需要说明的是，本文的模型尺度及实尺度的计算均采用壁面函数，因此，需要根据计算对象的雷诺数，对航行体及桨叶、桨毂表面第1层单元采用合适的高度，以使壁面y+值处于可适用壁面函数的范围，y+的具体数值见各算例。在阻力验证计算中，模型尺度与实尺度的计算网格拓扑关系以及各子域的单元类型均完全相同；而在敞水验证计算中，实尺度的计算网格由模型尺度网格按缩尺比（1∶1.5）放大得到。

1.2 计算设置

将航行体、推进器表面分别设为静止和旋转坐标系中的光滑无滑移壁面；计算域的下游边界面设为压力出口，其余边界面设为速度入口。动量方程和湍流输运方程采用二阶迎风格式离散，压力方程采用标准格式离散，压力与速度的耦合采用SIMPLE算法。

阻力计算采用定常模型；敞水及自航计算分别采用准定常模型和非定常（滑移网格）模型，在非定常计算的每个时间步前、后桨叶各转动1°。

实尺度自航计算在给定的航速下进行，通过调整转速（前、后桨相等），直至螺旋桨推力与航行体阻力之差小于后者的1%，即认为达到了航行体的自航点，结束计算。在自航的非定常计算中，推力、阻力均为一个桨叶旋转周期（360°）的时间平均值。敞水和自航的非定常计算收敛准则是，推力在连续2个旋转周期内的时间平均值之差小于0.5%。

2 计算精度验证

基于类似的航行体及对转桨的模型试验结果，分别验证阻力和敞水性能的数值计算精度；同时，进行实尺度计算，定性分析相关结果的合理性。因缺乏试验数据，对高速水下航行体的自航计算精度目前还无法验证。

2.1 航行体阻力

以图5所示的航行体为对象，验证阻力计算精度。在网格依赖性检验的基础上，先进行模型尺度的计算，对比试验数据以验证计算精度；然后，进行实尺度计算，并根据模型试验结果对实尺度阻力进行预报，间接验证实尺度阻力的RANS计算精度。

2.1.1 网格依赖性检验

为保证壁面y+值基本相同，采用相同的单元径向尺度，沿航行体轴向和周向按2∶2∶1的尺度比由疏到密建立3套网格M1，M2和M3，然后在模型尺度下进行了2个航速的阻力计算，雷诺数Rem分别为1.220×107和2.034×107，壁面y+的平均值均约为50。表1所示为3套网格参数的比较。航行体总阻力系数Ctm随网格尺度的变化如图6所示。随着网格的加密，Ctm趋于收敛；M2网格与M3网格Ctm的相对差仅为0.4%。因此，后续阻力计算采用M2网格的尺度，以兼顾计算精度和效率。

表1 3套阻力计算网格参数的比较Table 1 Comparison of the parameters for three sets of grids of resistance computation

2.1.2 模型和实尺度阻力计算精度验证

首先，验证模型尺度的阻力计算精度。表2所示为RANS计算结果与模型试验结果的比较。由于附体与壳体相比面积很小，因此分析形状因子时未区分壳体与附体，即

对模型试验，摩擦阻力系数Cfm按1957 ITTC公式计算：

对于RANS计算，Cfm直接采用计算结果。在计算雷诺数的范围内，模型总阻力计算值总体上低于试验值，但误差不超过3%；摩擦阻力的计算误差不超过3.5%，且随着航速的增高误差减小，这可能是因为计算是在全湍流状态下进行，而模型试验虽然采取了激流措施，但在航速较低时仍受到层流及过渡区的影响；形状因子的计算误差不超过4%。

实尺度网格拓扑与模型尺度相同，计算单元总数为1 985万。雷诺数Res=1.682×108，壁面y+平均值约为50。表3所示为实尺度航行体阻力的RANS计算结果与基于模型试验的预报结果间的比较，用于对实尺度阻力计算精度的间接验证。在预报中，忽略了形状因子的尺度效应，实尺度的形状因子(1+k)s为表2中模型试验雷诺数范围内(1+k)m的算术平均值，Cfs按1957 ITTC公式计算，Cts等于(1+k)s与Cfs的乘积，为航行体总阻力系数。在RANS计算结果分析中，Cfs和Cts均为直接计算值，(1+k)s为Cts与Cfs的比值。由表3可见，实尺度航行体总阻力的RANS计算值略低于基于模型试验的预报值，误差为3.3%。

表2 模型尺度航行体阻力RANS计算与试验的比较Table 2 Comparison of RANS-simulated model-scale resistance of the underwater vehicle with experimental data

表3 实尺度航行体阻力的RANS计算结果与基于模型试验的预报结果间的比较Table 3 Comparison between RANS-simulated full-scale resistance ofthe underwater vehicle and that predicted from model test results

2.2 对转桨敞水性能

以某高速航行体对转桨为对象，验证敞水性能的RANS计算精度。该对转桨前桨11叶、后桨9叶，前、后桨的直径之比为1.062。首先，进行模型尺度的网格依赖性分析，确定合适的网格尺度；然后，进行模型尺度和实尺度敞水性能的计算与试验比较。

2.2.1 网格依赖性检验

计算在模型尺度进行，进速系数J=1.0和1.4时，如式（3）所示的0.7R（R为桨直径）处雷诺数分别为8.45×105和9.05×105。

式中：VA为进速；b0.7R，n和D分别为桨叶0.7R处的弦长、转速与直径；ν为水的运动粘性系数。由疏到密建立了3套网格G1，G2和G3，其单元尺度之比为2∶2∶1，面单元和体单元的密度分布及增长率保持不变。表4所示为3套网格的参数比较，其中面单元相对尺度为面单元尺度与前桨直径的比值。

表5所示为网格尺度不同引起推力系数KT和扭矩系数KQ计算结果的相对变化。表中KT和KQ均以前桨的直径无量纲化，下标F和A分别表示前、后桨。分别用R1，R2和R3代表G1，G2，G3的KT或KQ计算结果，则表中C21=R2/R1-1，表示R2相对R1的变化，C32=R3/R2-1，表示R3相对R2的变化。从表中可以看出，前、后桨的推力和扭矩基本上是随网格尺度的减小而减小；3套网格计算结果之差均在1.2%以内。经权衡精度与效率，后续的计算采用G2网格的尺度。

表4 3套敞水计算网格参数的比较Table 4 Comparison of the parameters for three sets of grids of open-water computation

表5 敞水性能计算结果随网格尺度的相对变化Table 5 Influence of grid size on computed open-water performance

2.2.2 敞水性能计算

用于验证的对转桨模型与实桨的尺度比为1∶1.5，实桨的计算网格由桨模网格放大得到。在计算工况范围内，桨模与实桨在0.7R处的雷诺数分别为8.45×105～9.05×105和3.17×106～3.39×106，桨叶表面的y+平均值分别为35～50和130～200，均在壁面函数的适用范围内。当前、后桨的叶数相互不构成整倍数关系时，两者相互作用的非定常性很弱［14］，非定常计算结果的时间平均值与准定常结果很接近，因此采用准定常计算以节约计算时间。

图7所示为敞水性能的模型和实尺度RANS计算结果与模型试验结果比较，图中η0为敞水效率，下标F和A分别表示前、后桨。在计算进速系数范围内，总推力系数和总扭矩系数的计算误差分别小于2%和4%，表明模型尺度RANS计算具有较高的精度。比较实尺度与模型尺度的RANS计算结果发现，KT的变化很小，而KQ则有较明显的降低，导致实桨的敞水效率比模型桨的高2.5%左右，该结果与单桨敞水性能尺度效应定性一致，表明对转桨实尺度敞水性能的RANS计算方法比较可靠。

3 实尺度自航计算与分析

基于上述计算方法，对采用对转桨推进的某高速水下航行体进行实尺度阻力、敞水及自航计算，并分析其推进因子。对转桨的前桨9叶、右旋，后桨7叶、左旋，二者转速相同。

3.1 航行体阻力计算

阻力计算单元总数约1 127万。计算航速范围为设计航速的64%～109%，雷诺数为8.07×107～13.83×107，航行体壁面y+的平均值约为50。图8所示为该航行体实尺度阻力的RANS计算模型，图9所示为拖航阻力系数计算结果。图10所示为航行体摩擦阻力系数的RANS结果与1957 ITTC公式计算结果的比较，由图可见，前者略高，两者间的差异随航速的增大而减小，在设计航速下相差1.9%。图11所示为形状因子(1+k)s随航速的变化，虽然形状因子随航速的增高而变大，但在计算航速范围内其差值小于1%。

3.2 对转桨敞水性能计算

敞水计算单元总数约为4 634万，桨叶0.7R处雷诺数范围为2.38×106～2.57×106，桨叶表面y+的平均值约为210。分别采用准定常和非定常方法模拟前、后桨之间的相互作用，其中准定常方法用于模拟桨叶旋转过程中某一瞬间的流动，而非定常方法则是在时域中模拟桨叶旋转过程中的流动。图12和图13所示分别为实尺度对转桨敞水性能的RANS计算模型和计算结果，其中非定常方法模拟的结果为桨叶旋转一周的平均值，图中扭矩系数KQ的第2下标0表示敞水，以与后文中的自航数据相区别。从图13中可以看出，采用准定常方法与非定常方法模拟的敞水性能结果非常接近，差值在1%以内。

3.3 自航计算结果

自航计算在设计航速下进行，计算单元总数为4 827万，航行体雷诺数为12.69×107，壁面y+的平均值约为50。分别采用准定常和非定常方法进行计算，自航转速根据1.2节所述方法确定。图14所示为自航的RANS计算模型。为了节约计算时间，非定常计算采用准定常计算得到的自航转速及流场作为初始值。图15所示为非定常计算的轴向力收敛历程。图中，KQBF和KQBA分别为前、后桨的扭矩系数。其中推力系数和扭矩系数的脉动幅值分别为各自平均值的0.2%与0.1%，而Cts的脉动幅值则约为平均值的2.5%。Cts的脉动主要源于由桨叶尾流引起的航行体尾端面压力脉动。

自航的非定常与准定常计算结果间的比较见表6，表中Ns为航行体推力和阻力平衡时前、后桨的转速（两者相等）。由于航行体伴流不均匀性的影响，自航非定常计算与准定常计算结果间的差别比敞水情况下的大，主要表现为非定常计算的后桨负荷比准定常的高。从理论上讲，非定常计算模型更接近于物理实际，因此在前、后桨扭矩平衡要求较高的设计场合，采用非定常计算相对可靠。

表6 自航计算的准定常与非定常计算结果的相对差Table 6 Relative difference between quasi-steady and unsteady results of self-propulsion simulation

图16所示为自航状态基于Q准则的尾部涡流形态（准定常结果）。从图中可以看出，在桨叶边缘、附体的尾缘和航行体尾端面均出现了明显的涡流，较合理地反映了桨叶对周围流场的扰动。由于前、后桨旋转域间交界面上的网格无法做到相同，因此前桨尾涡在进入后桨域后因网格耗散几乎消失了。

图17所示为在尾端面下游0.05倍前桨直径处，横截面内切向速度沿半径的分布。其中，切向速度为沿圆周线的平均值，并用当地半径桨叶的旋转线速度进行无量纲化。所谓当地半径，是指在航行体中纵剖面上，将通过前、后桨盘半径外端点的线段延长到该纵向位置，线段终点与航行体轴线的距离。图17中的横坐标为径向相对坐标2r/DF，其中r为径向坐标值。0.65R至叶梢的切向速度几乎为零，但内半径的切向速度仍有残留。总体来看，在设计航速下，该航行体尾流基本没有旋转，表明前、后桨的扭矩平衡较好。

3.4 自航因子分析

综合阻力、敞水及自航计算结果，采用等推力法进行自航因子分析，结果如表7所示。表中：J0为根据KT在图13对应的推力系数曲线上插值得到的进速系数；η0为J0所对应的效率曲线上的值；KQ的第2下标B表示自航状态，用于与敞水状态的0相区别；ηR为相对旋转效率；ηD为推进效率；下标F和A分别表示前桨和后桨。

表7 自航分析结果Table 7 Results of self-propulsion analysis

分析比较表明：

1）自航因子计算值处于合理的范围。参考文献［15］，水下航行体的伴流分数w和推力减额系数（1-t）的范围分别为0.10～0.25和0.82～0.90，本文中航行体的计算结果也处于此范围内。

2）前桨的伴流分数高、后桨的伴流分数低，分析原因认为主要为：一是前桨更靠近尾附体，附体的尾流速度亏缺量更大；二是前、后桨毂径之比为1.31，而直径之比仅为1.06，这就意味着前桨盘面有较多的面积处于桨毂边界层中。此外，桨叶的抽吸作用会导致桨盘面处的实效伴流低于标称伴流，如果这种作用在桨盘前较弱、桨盘后较强，也会成为原因之一，不过该问题还有待进一步的研究。

3）推进因子及推进效率的准定常与非定常计算结果间的误差在2%以内。若采用准定常计算方法得到的推力减额较低、推进效率较高，则有可能导致设计桨偏重，预报的航速偏高。

4 结 论

本文通过求解RANS方程，对采用对转桨推进的高速水下航行体开展了实尺度自航计算与分析研究。通过与模型试验结果的比较，验证了在航行体阻力和对转桨敞水性能方面，采用本文方法计算精度良好；间接的验证及定性分析表明，实尺度阻力和敞水的计算结果合理、计算方法可靠，但实尺度自航计算仍有待验证。对某航行体实尺度的自航计算与分析表明：

1）采用准定常和非定常方法计算得到的推进因子其数值比较接近，也在合理范围内，可望为高速水下航行体对转桨设计提供较可靠的输入，从而提高设计精度、减少模型试验，缩短设计周期。

2）在自航模拟中，准定常计算方法具有精度合理、计算量小的优点，比较适合工程应用；但当航速预报精度以及前、后桨的扭矩平衡精度要求较高时，采用非定常计算方法更合理。