0 引　言

对于工作在船后或艇后的螺旋桨，由于其处在径向和周向都不均匀的三维伴流场中，因而会产生周期性变化的非定常力，进而激发出离散谱噪声。所以，螺旋桨非定常力的分析和研究对预报螺旋桨线谱噪声和提高潜艇的隐身性至关重要。对螺旋桨非定常力的研究，模型试验是比较传统的方法，但由于试验周期长、成本高，还会受到尺度效应、外界干扰、伴流场的模拟等因素的影响。因而模型试验方法受到了很大的限制，这方面的文章也很少见。从20世纪70年代开始，基于势流理论的螺旋桨非定常力数值分析得到了迅速发展。在这方面的代表学者，国外主要有Kerwin^[1]，Hoshino^[2] ，国内有陈家栋^[3]、谭延寿^[4]、熊鹰^[5]等。但这种方法需要处理库塔条件和假定尾涡形状，并且很难较好地处理粘性的影响。

随着计算流体力学CFD技术的迅猛发展，粘性数值计算方法逐渐成为研究螺旋桨非定常力的另一有效方法。张志荣等^[6]基于CFD技术，采用3种方法对斜流条件下大侧斜螺旋桨的非定常水动力性能进行预报研究，分析了不同方法对计算螺旋桨非定常水动力性能的适用性。胡小菲^[7]对非均匀来流下的螺旋桨所受的非定常力进行数值模拟，并对网格的收敛性进行分析，计算结果与试验值吻合较好。黄振宇^[8]采用SST k-ω湍流模型和滑移网格技术对非均匀来流条件下的2个螺旋桨所受的非定常力进行数值模拟，计算结果与试验值吻合较好。在此基础上，对全附体潜艇模型后的大测斜螺旋桨的非定常力进行计算。刘登成等^[9]基于Fluent软件，以DTMB4119为研究对象，对速度进口到桨盘面的距离以及网格数量进行分析。并根据建立的数值方法，对HSP螺旋桨非定常力进行数值预报。马超^[10]采用滑移网格技术，结合RNG湍流模型对Suboff全附体艇体下的螺旋桨的非定常力进行数值计算，获得螺旋桨转动一周过程中轴承力的变化规律。

本文基于STAR-CCM+仿真软件，使用滑移网格技术，对非均匀来流下的螺旋桨非定常力预报方法进行研究，通过与试验值进行比较，具体分析了时间离散格式、螺旋桨每时间步旋转的度数以及网格密度对计算结果的影响。并在此数值方法的基础上，对艇体非均匀伴流场下的螺旋桨所受到的非定常力进行了数值分析。

1 控制方程与湍流模型的选取 1.1 控制方程

不可压缩粘性流体的基本控制方程由连续性方程和RANS方程组成，其张量形式为：

$\frac{{\partial {u_i}}}{{\partial {x_i}}} = 0\text{，}$

$\begin{align}&\frac{\partial }{{\partial t}}(\rho {u_i}) + \frac{\partial }{{\partial {x_j}}}(\rho {u_i}{u_j}) = - \frac{{\partial p}}{{\partial x{}_i}} + \frac{\partial }{{\partial {x_j}}}\\&\left[ {{\mu _0}(\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}) - \frac{2}{3}{\mu _0}\frac{{\partial {u_l}}}{{\partial {x_l}}}{\delta _{ij}}} \right] + \frac{\partial }{{\partial {x_j}}}( - \rho \overline {u_i'u_j'} ) + \rho {f_i}\text{。}\end{align}$

式中：ρ为流体密度；p为静压；f_i为单位质量的质量力；u_i，u_j，u_l均为速度分量；μ₀为分子粘性系数； $ - \rho \overline {u_i'u_j'} $ 为由于时均化而增加的湍流雷诺应力项。

1.2 湍流模型

显然，上面的控制方程不封闭，因而要对流场进行求解还需要选取适当的湍流模型。通过对相关文献的阅读，发现对于螺旋桨非定常水动力性能的数值模拟，研究人员较多的采用SST k-ω湍流模型。因此，本文选取SST k-ω湍流模型对方程组中的雷诺应力项进行封闭。

该模型由Menter提出，又称为剪切应力输运k-ω模型，它综合了边界层外部k-ε模型独立性和近壁k-ω模型稳定性的优点。在近壁面区域采用k-ω湍流模式，在边界层外部采用k-ε湍流模式。该模型的k方程和ω方程分别为：

$\frac{\partial }{{\partial t}}(\rho k) + \frac{\partial }{{\partial {x_i}}}(\rho k{u_i}) = \frac{\partial }{{\partial {x_j}}}\left[ {(\mu + \frac{{{\mu _t}}}{{{\sigma _k}}})\frac{{\partial k}}{{\partial {x_j}}})} \right] + {P_k} - {\beta '}\rho k\omega \text{，}$

$\frac{\partial }{{\partial t}}(\rho \omega ) + \frac{\partial }{{\partial {x_i}}}(\rho \omega {u_i}) = \frac{\partial }{{\partial {x_j}}}\left[ {(\mu + \frac{{{\mu _t}}}{{{\sigma _\omega }}})\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + \alpha \frac{\omega }{k}{P_k}{\rm{ - }}\beta \rho {\omega ^2}\text{。}$

2 非均匀来流下螺旋桨非定常力的数值模拟 2.1 计算对象与计算域

对于非均匀来流下螺旋桨非定常力的计算，本文选取DTMB4119桨为研究对象。该桨模为3叶桨，直径D=0.305 m，毂径比D_h/D=0.2，设计进速系数J=0.833，无侧斜和纵倾，右旋桨，几何模型如图1所示。

在螺旋桨上建立柱坐标系（x，r，θ），原点为桨盘面与桨轴的交点，x轴与桨轴重合，指向下游为正；r为径向，指向外为正；θ为周向，从桨后向桨前看顺时针为正，0°在12点钟位置。

考虑到要使用滑移网格技术进行非定常计算，因而计算流场区域分为两部分，外部的静止域和内部的旋转域。外部区域取为桨前4D、桨后7.5D、半径为2.5D的大圆柱形，内部区域取为桨前0.23D、桨后0.4D、半径为0.65D的小圆柱形，2个圆柱的轴线与桨轴重合，整个计算域如图2所示。

2.2 边界条件与网格划分

整个计算区域的边界分为进流边界、出流边界、远场边界、螺旋桨以及桨轴表面。具体的边界设置为：进流边界取为速度入口边界类型，出流边界取为压力出口边界类型，远场边界取为对称面边界类型，螺旋桨以及桨轴表面取为无滑移壁面边界类型。对于非均匀来流的设置，本文采用Jessup（1990）^[11]在DTRC报告中九阶角频的轴流场，在进流边界各个半径处（R为螺旋桨半径）的周向速度分布如图3所示。本文在设计工况J=0.833下，取螺旋桨转速为25rps，则平均来流速度U₀=6.35 m/s。

STAR-CCM+自带的网格生成技术可以快速高效的生成质量很好的计算网格，本文采用以六面体为核心的切割体网格（Trimmer）、拉伸型网格（Extruder）和棱柱层网格（Prism Layer Mesher）对计算域进行离散。考虑到在入口边界处要进行非均匀来流的设置，为了更好的反应各个半径周向不同位置处速度的变化，对入口边界的网格划分，要从外到内逐步加密，如图4所示。最终得到的入口边界速度分布云图如图5所示，从图中可以看出，速度分布的九阶特征比较明显，因而入口边界的网格满足非均匀来流速度设置的要求。另外，根据螺旋桨周围流场的流动特点，通过体积控制（Volumetric Controls）功能模块对靠近螺旋桨的区域分层逐步加密（见图6），以便更好地捕捉流动细节，获得更加精确地流场信息。最后，考虑到在螺旋桨的导边、随边、叶梢以及叶根处各个物理量的变化梯度较大，对其也进行了加密（见图7）。

2.3 计算结果与分析

基于上述的计算模型、网格系统，运用有限体积法对螺旋桨的非定常力进行数值模拟。通过将数值计算结果与Jessup （1990）^[1]的试验数据进行对比，分析时间离散格式与单位时间步长旋转角度以及网格密度对计算结果的影响，得到了适用于螺旋桨非定常力计算的数值分析方法。

2.3.1 时间离散格式与单位时间步长螺旋桨旋转角度的确定

为了研究时间离散格式与单位时间步长螺旋桨旋转角度对非定常力的影响，本文分别在1阶时间离散格式和2阶时间离散格式下，保持网格划分方式以及网格密度不变，通过改变螺旋桨单位时间步长的旋转角度（3.6°，2.7°，1.8°，0.9°，0.45°）对螺旋桨的非定常力进行数值计算。将得到的数值计算结果与试验值进行比较，如图8所示。图中，0°在12点钟位置，从桨后向桨前看，顺时针方向旋转。

通过比较分析，可以得出如下结论：首先，无论是采用1阶时间离散格式还是2阶时间离散格式，随着螺旋桨单位时间步长旋转角度的减小，数值计算结果与试验值越来越接近，精度也越来越高。另外，在螺旋桨单位时间步长旋转角度相同的情况下，2阶格式的数值计算结果与试验值的吻合度明显比1阶格式更高；最后，对于2阶时间离散格式，当螺旋桨单位时间步长的旋转角度小于1.8°时，数值计算结果与1.8°差别很小。综上所述，本文认为对于螺旋桨非定常力的计算，螺旋桨单位时间步长的旋转角度应为1.8°同时选取2阶时间离散格式。

2.3.2 网格密度对螺旋桨非定常力的影响

对于非均匀来流下螺旋桨非定常力的数值模拟，由于靠近桨叶表面附近的流动较为复杂，各个物理量的变化梯度较大，因而要对其周围的网格进行加密。而在STAR-CCM+仿真软件中，可以通过设置桨叶表面网格的边界增长率（Boundary Growth Rate）的快慢，对其周围网格疏密程度进行调整。为了研究桨叶表面周围网格疏密程度对计算结果的影响，本文设计了3套计算网格（见图9），边界增长率分别为medium，slow，very slow，对应的旋转域的网格数为90万、130万和250万。根据2.3.1节的研究结论，在2阶时间离散格式以及螺旋桨单位时间步长的旋转角度为1.8°的情况下，对3套计算网格下的螺旋桨非定常力进行数值模拟，得到的计算结果如图10所示。

从图中可以看出，在给定的螺旋桨单位时间步长的旋转角度和时间离散格式下，当桨叶表面网格的边界增长率设定为medium时，得到的非定常力数值计算结果与试验值吻合的较好；并且继续对桨叶表面周围的网格进行加密，得到的非定常力数值计算结果与medium的结果相比差别很小。因而对于螺旋桨非定常力的计算，桨叶表面网格的边界增长率设定为medium就可以达到计算精度的要求。

3 全附体潜艇模型后的螺旋桨非定常力计算

基于上节的数值方法和旋转域网格划分形式，本节对Suboff全附体潜艇模型下3个螺旋桨（DTMB4119，DTMB4381，DTMB4383）的非定常力进行数值模拟。Suboff全附体潜艇模型由一个水滴形轴对称体、一个指挥台围壳和4个尾翼组合而成，具体的几何参数见文献[12]。在对自航状态进行计算之前，本节先对无桨时的全附体潜艇尾流场进行了数值分析，得到了桨盘面r/R=0.25（R为主艇体最大半径）处的速度分布，并与试验值进行比较（见图11）。从图中可以看出，在幅值方面，各个方向上的伴流分数与试验值吻合的较好；对于相位，从0°到180°，计算值与试验值的相位吻合较好，而从180°到360°，数值模拟值与试验值的相位有一定角度的偏差。由于数值模拟计算结果的对称性较好，x方向伴流分数的峰值都与艇后舵的位置相对应，且最大值对应于指挥台围壳的位置，因而可以认为这二者的偏差是由试验误差造成的。总体而言，计算值与试验值吻合的较好，验证了潜艇网格划分的合理性，为艇桨结合计算打下了基础。

4381和4383均为5叶桨，毂径比均为0.2，但前者无侧斜，后者的侧斜角为72°。3部螺旋桨的直径统一取为 $D = 0.4{D_H}$ ，D_H为Suboff主艇体最大直径，螺旋桨转速统一设置为720 r/min。将计算得到的单个桨叶非定常力变化规律与刘志华^[13]的计算结果进行比较如图12所示。从图中可以看出，在相位方面，本文计算值与其计算结果吻合较好，非定常力的变化规律与桨盘面处轴向伴流分数的分布相对应，伴流分数的低点对应非定常力的高点，伴流分数的高点对应非定常力的低点。在幅值方面，本文计算的非定常力变化幅度与刘的计算结果相比稍微有点偏大，在最低点和最高点处这种差别更明显。因而，对于艇桨结合下螺旋桨的非定常力数值计算，迫切需要相应的试验结果进行验证。另外，从图中还可以看出，在计算条件和网格划分都相同的条件下，大侧斜螺旋桨可以明显降低单个桨叶非定常力的变化幅度。

图13给出了3部螺旋桨无因次化的全桨非定常力随旋转角度的变化情况，从图中可以直观看出，推力的变化幅度随着桨叶叶数的增加以及侧斜的增大而降低。为了对非定常力的变化进行定量研究，本文将计算得到的全桨非定常力进行傅里叶分析，得到了1阶叶频的轴向轴承力，并与拟定常方法的结果进行比较如图14所示。对于CFD方法，通过傅里叶分析得到的（4119，4381，4383）3部螺旋桨的轴向轴承力的1阶叶频幅值分别为2.267 N，1.56 N，0.467 N。其中，5叶桨4381的幅值相比3叶桨4119降低了30% ，大侧斜螺旋桨4383的幅值相比无侧斜桨4381降低了70%。而采用拟定常经验公式得到的不同螺旋桨的轴向轴承力变化不大，基本相当，因而拟定常方法不适应多桨叶螺旋桨尤其是大侧斜螺旋桨的轴承力计算，无法满足工程需要。

4 结　语

本文基于STAR-CCM+对非均匀来流下的螺旋桨非定常力进行了数值模拟，研究了时间离散格式、单位时间步长螺旋桨旋转角度与网格密度对计算结果的影响。并在此基础上，对Suboff全附体潜艇模型下3个螺旋桨的非定常力进行数值模拟。得到结论主要有：

1）时间离散格式和单位时间步长螺旋桨旋转角度对非均匀来流下螺旋桨非定常力的影响较大，通过与试验结果的比较分析，认为螺旋桨单位时间步长的旋转角度应为1.8°同时选取2阶时间离散格式。

2）在网格密度方面，桨叶表面网格的边界增长率设定为medium就可以达到计算精度的要求。

3）对于全附体潜艇模型下螺旋桨的非定常力计算，桨叶侧斜以及桨叶个数对非定常力的变化幅度影响较大，随着侧斜角和桨叶个数的增加，非定常力的变化幅度会大幅降低。