斜流中艇后螺旋桨水动力数值计算方法

张文照，肖昌润

(1.海军工程大学,湖北 武汉 430033； 2.海军后勤技术装备研究所,北京 100072)

0 引 言

斜流条件下，艇后螺旋桨水动力性能研究一直是进行潜艇操纵性研究的关键，目前主要通过模型试验的方法来进行研究。但模型试验存在尺度、伴流场的模拟、结果换算和经费等众多问题的困扰，受到极大限制。数值方法可以克服上述试验方法存在的不足，势流方法是螺旋桨性能分析中常用的数值方法。在非均匀流中螺旋桨性能分析也常采用势流方法，比较突出的有Kinnas[1]和Hosino[2]。国内有陈家栋[3]、熊鹰[4]、刘小龙[5]等用势流方法分析给定不均匀流场中螺旋桨性能。

但采用势流方法处理艇后螺旋桨这种多体结构的水动力问题变得十分困难，并且势流方法在处理非均匀进流中螺旋桨性能时，需要处理库塔条件和假定尾涡形状，有一定的人为因素。而粘流CFD方法则更方便处理复杂的几何模型，实现对流场的真实模拟。Funeno[6]用粘流CFD方法分析了大侧斜桨非定常性能，WATANABE[7]使用商业软件Fluent研究了非均匀伴流中的螺旋桨性能。当前粘性CFD的迅速发展，使得用CFD方法研究艇后螺旋桨的水动力性能变得越来越容易。本文应用Fluent软件计算了斜流条件下艇后螺旋桨的水动力，为开展潜艇操纵性研究奠定了基础。

1 研究对象

本文所研究的潜艇模型为加拿大DRDC潜艇水动力研究计划中的通用艇模[8]，模型长4.36 m，直径0.49 m。螺旋桨模型为5叶桨，直径0.22 m，盘面比0.5。图1给出了艇体带桨的几何模型。

图1 计算所用带桨潜艇模型Fig.1 Submarine model

2 计算方法描述

本文采用基于粘流理论的数值方法模拟艇后螺旋桨周围的流动，以商用软件Fluent为求解器求解不可压三维RANS方程计算艇后螺旋桨的推力和扭矩。

2.1 控制方程

螺旋桨水动力计算所选取的坐标系为固定在螺旋桨上的旋转直角坐标系，定义坐标轴x为螺旋桨的旋转轴，指向下游为正。控制方程采用Reynolds平均的NS(RANS)方程，其不可压缩的无量纲形式为：

(1)

式中：ui和uj为绝对速度各分量；p为压力；Re为雷诺数；

(4)

式中：νt为湍流运动引起的涡运动粘性系数；

采用SSTk-ω湍流模型来计算雷诺应力项，以实现控制方程组的封闭。SSTk-ω湍流模型混合了k-ω模型计算近壁面区域粘性流动的可靠性和k-ε模型计算远场自由流动的精确性，并在螺旋桨性能计算中反复验证其较高的精度。

控制方程采用有限体积法离散，有限体积法具有良好的守恒性，能够采用各种形式的网格以适应复杂的几何轮廓。控制方程中的对流项采用二阶迎风格式离散，速度与压力的耦合迭代采用SIMPLE 算法。离散后形成的代数方程用逐点Gauss-Seidel迭代法求解，并且采用代数多重网格方法加快求解的收敛速度。

2.2 滑移网格法

在进行艇后螺旋桨水动力计算时，螺旋桨旋转，艇体(包含附体)静止。对于这种问题，Fluent软件有混合面法、多参考系法、滑移网格法和动网格法4种处理方法，前2种方法都假设流动为定常，后2种方法都假定流动为非定常。动网格法需要配合诸如层生成吞并等技术，设置复杂，而滑移网格法在处理定转子问题上更易于实现。本文对于艇后螺旋桨周围流场的计算采用滑移网格法，包含艇体和附体的部分设置为静止区域，螺旋桨部分为滑动区域，计算网格随螺旋桨转动。在2个区域的交界面上物理量的传递采用插值方法处理。

2.3 计算域与网格划分

计算域的进口边界设置为艇首上游0.8L处(L为艇模长度)的半球面，出口边界设置在螺旋桨下游2L处，柱面外边界也为进口，柱面半径取1L。整个计算域分割成两大部分，如图2所示，包含螺旋桨的那一部分区域为旋转区域，其他区域为静止区域。

图2 计算区域示意图Fig.2 Computing field

包含螺旋桨的旋转部分由于螺旋桨复杂的几何外形，划分为六面体网格较难实现，所以对旋转区域仍采用四面体网格划分(见图3)。

静止和旋转区域交界面上的网格划分见图4。

图3 旋转区域的网格划分Fig.3 Mesh grid of rotation field

图4 静止和旋转区域交界面处的网格划分Fig.4 Mesh grid of the interface of static field with rotation field

包含吊艇体、围壳和十字舵的静止区域在进行合理的分块后较易实现六面体网格划分。静止区域的分块如图5所示，分块会给网格划分带来一个问题，即指挥台围壳和尾翼表面一层的网格尺寸会延续到整个计算域的外边界，然而在计算域的外边界，其他方向上的网格划分已经很疏，这样就会使得网格长宽比过大。于是，在网格划分用如图6所示的处理方法，即采用“斜向扩展”的方式来防止网格在计算域外边界上长宽比过大。

图5 静止区域网格划分Fig.5 Mesh grid of static field

图6 指挥台围壳网格斜向扩展法Fig.6 Oblique extension method

图7 围壳和尾翼端面网格划分Fig.7 Mesh grid of shell and end face

此外，对于围壳和尾翼端面网格的划分，将上端面先分成如图7所示的2个端部面和1个中间面，当支柱上端面边的节点固定后，这样的分块能有效控制网格层数。

2.4 边界条件和初始条件的设置

计算域的进口为图2中右侧的圆面和径向的柱面，进口设置为速度进口，给定来流速度的大小和方向。计算域的出口为左侧的圆面，出口设置为压力出口，给定压力大小。螺旋桨、艇体、指挥台围壳和尾附体表面设为物面无滑移条件，静止区域与旋转区域相重叠的面设置为交界面，完成2个区域的数据传递。用进口速度初始化整个流场。

3 计算结果的比较与分析3.1 0°攻角不同进速非定常计算结果

首先分别计算了艇体0°攻角下，螺旋桨艇后进速Jb为0.8，1.0，1.2，1.4，1.6工况下，螺旋桨转动1周过程中推力和扭矩的变化值。选取设计点Jb=1.0，画出该工况下螺旋桨转动1周所产生的推力和扭矩的变化曲线(见图8和图9)。通过曲线可以看出螺旋桨转动1周，推力和扭矩呈现出规律的周期性，周期数等于螺旋桨的页数。螺旋桨转速通常很高，因此推力和扭矩的这种周期性的脉动相对于潜艇运动而言十分快速，因此在研究潜艇操纵运动时，给定艇体运动状态的情况下(如潜艇航速、纵倾角等)，可以忽略螺旋桨推力和扭矩的这种周期脉动，而用一个周期的平均值作为艇体该运动状态下螺旋桨产生的推力和力矩计算值。

表1和表2分别列出了螺旋桨艇后进速Jb为0.8，1.0，1.2，1.4，1.6工况下计算得到平均推力系数Kt和平均扭矩系数10Kq，并和试验值做比较，计算误差在3%以内。

图8 Jb=1.0推力系数1周内的变化Fig.8 Trust coefficient in a period when Jb=1.0

图9 Jb=1.0扭矩系数1周内的变化Fig.9 Torque coefficient in a period when Jb=1.0

进速JbKt计算值Kt试验值相对误差/%080275028019100242024408120211020816140175017123160136013325

表2 不同进速螺旋桨平均扭矩系数

3.2 艇体攻角对螺旋桨水动力性能影响计算

为了研究艇体攻角对艇后螺旋桨推力和扭矩的影响，计算设计进速Jb=1.0， 艇速u=0.7 m/s，艇体处于0°，5°，10°，15°，30°攻角工况下艇后螺旋桨平均推力、扭矩系数(见图10)。由于整个计算对象几何对称特性，负攻角与正攻角下的桨模水动力对称。

与试验值相对比，计算值的相对误差随着攻角的增大而增大，总体控制在7%以内，该误差是由于在大攻角时，艇体流体分离比较严重，而计算网格不能有效地捕捉到该现象导致的。

图10 Jb=1.0时艇后螺旋桨推力、扭矩系数随艇体攻角的变化Fig.10 Behind-boat trust of the propeller with different pitch angle when Jb=1.0

从图10可看出，螺旋桨推力和扭矩是随着艇体攻角呈现出先减小后增大的趋势。究其原因，是因为当艇体产生攻角时，螺旋桨进流速度会减小，伴流分数会增大，当攻角比较小的时候伴流分数的增大导致的螺旋桨实效进速的提高占主要作用，因此在小攻角时螺旋桨推力和扭矩均随着攻角的增加而降低。但是，随着攻角的不断变大，攻角导致的螺旋桨进流速度的降低越来越严重，进而超过伴流分数减小的影响，此后，螺旋桨的推力和扭矩随着攻角的增大而有所增加。

4 结 语

本文采用CFD方法对艇后螺旋桨的水动力性能进行计算，并对计算结果进行比较和分析后得出结论如下：

1)本文采用的数值计算方法对带攻角的艇后螺旋桨推力和扭矩是有效的；

2)随着攻角的增大，螺旋桨推力和扭矩的计算误差有所增加，总体控制在7%以内；

3)在同一进速下，艇后螺旋桨推力和扭矩随着艇体攻角的增大呈现出先减小后增大的趋势。

[1] KINNAS S A,HSIN C J, KEENAN D.A potential based panel method for the unsteady flow around open and ducted propellers[C].Proceedings of the 8thSymposium on Naval Hydrodynamics,1990.

[2] HOSHINO T.Hydrodynamic analysis of propeller in unsteady flow using a surface panel method[C].Journal of Society of Naval Architectures of Japan,1993.

[3] 陈家栋,董世汤.非定常螺旋桨表面压力面元法预报计算[J].中国造船，1998(1)：21-26.

[4] 熊鹰.非均匀流中螺旋桨空泡和脉动压力地数值和试验研究[D]，武汉：武汉理工大学，2002.

[5] 刘小龙,王国强.螺旋桨非定常性能预估的面元法[J].船舶力学,2006,10(2).

[6] FUNENO I. Analysis of unsteady viscous flows around a highly skewed propeller(in Japanese)[J]. Kansai Society of Naval Architects,2002,237：39-45.

[7] WATANABE T.Simulation of steady and unsteady cavitation on a marine propeller using a RANS code，Fifth International[J].Symposium on Cavitation (CAV2003)，Osaka, Japan,2003：1-4.

[8] FOURNIER E.Wind tunnel investigation of the DREA po-wered MkI submarine model[R].CANADA:DREA,1993.