不均匀进流场下螺旋桨非定常力的数值模拟

张凯奇，叶金铭，于安斌

(海军工程大学 舰船工程系，湖北 武汉 430033)

0 引 言

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

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

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

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

1.1 控制方程

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

1.2 湍流模型

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

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

2 非均匀来流下螺旋桨非定常力的数值模拟

2.1 计算对象与计算域

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

图 1 DTMB4119桨的几何模型Fig. 1 DTMB4119 model

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

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

图 2 数值模拟计算域Fig. 2 Computation domain of numerical simulation

2.2 边界条件与网格划分

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

图 3 DTRC 4119桨的进流条件Fig. 3 The inflow conditions of DTRC 4119

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

图 4 入口边界网格划分Fig. 4 Inlet boundary mesh

图 5 入口边界的速度云图Fig. 5 velocity contour of inlet boundary

图 6 对称面上的网格分布Fig. 6 Mesh distribution of symmetry plane

2.3 计算结果与分析

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

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

图 7 桨叶和桨榖表面上的网格分布Fig. 7 Mesh distribution of blades and hub

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

图 8 时间离散格式及旋转角度对非定常力的影响Fig. 8 The impact of temporal discretization and rotation angle on unsteady hydrodynamics

通过比较分析，可以得出如下结论：首先，无论是采用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就可以达到计算精度的要求。

图 9 三套计算网格Fig. 9 3 sets of mesh

图 10 网格疏密程度对非定常力的影响Fig. 10 The impact of mesh on unsteady hydrodynamics

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

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

图 11 桨盘面处的速度分布Fig. 11 Velocity distribution at propeller disk

4381和4383均为5叶桨，毂径比均为0.2，但前者无侧斜，后者的侧斜角为72°。3部螺旋桨的直径统一取为，DH为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%。而采用拟定常经验公式得到的不同螺旋桨的轴向轴承力变化不大，基本相当，因而拟定常方法不适应多桨叶螺旋桨尤其是大侧斜螺旋桨的轴承力计算，无法满足工程需要。

图 12 单个桨叶的非定常力Fig. 12 Unsteady hydrodynamics of single blade

图 13 全桨的非定常力Fig. 13 Unsteady hydrodynamics of blades

图 14 1阶叶频轴向轴承力Fig. 14 Axial bearing force of 1-st order blade frequency

4 结 语

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

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

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

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

[ 1 ]KERWIN J E, LEE C S. Prediction of steady and unsteady marine propeller performance by numerical lifting-surface theory[J]. Societyof Naval Architects and Marine Engineers Transactions, 1978, 86: 218–253.

[ 2 ]HOSHINO T. Hydrodynamic analysis of propeller in unsteady flow using a surface panel method[J]. Journal of the Society of Naval Architects of Japan, 1993, 74: 71–87.

[ 3 ]陈家栋, 董世汤. 非定常螺旋桨表面压力面元法预报计算[J].中国造船, 1998(1): 21–26.CHEN Jia-dong, DONG Shi-tang. Prediction of pressure distributions on unsteady propeller blades using potential-based panel method[J]. Ship Buiding of China, 1998(1): 21–26.

[ 4 ]谭延寿, 贺伟. 螺旋桨非定常轴承力计算[J]. 船海工程, 2006,171(2): 42–45.TAN Yan-shou, HE Wei. Calculation of unsteady shaft forces of propeller[J]. Ship & Ocean Engineering, 2006, 171(2):42–45.

[ 5 ]熊鹰. 非均匀流中螺旋桨空泡和脉动压力的数值和试验研究[D]. 武汉: 武汉理工大学, 2002.XIONG Ying. Numerical and experimental research on propeller-induced pressure fluctuations and cavitation in nonuniform flow[D]. Wuhan: Wuhan University of Technology,2002.

[ 6 ]张志荣, 洪方文. 斜流中螺旋桨性能数值分析方法研究[C]//船舶力学学术会议暨《船舶力学》创刊十周年纪念学术会议. 银川, 2007.ZHANG Zhi-rong, HONG Fang-wen. Research on performance analysis method for propeller in oblique flow[C]//Proceedings of Ship Dynamics Meeting, 2007, 9.

[ 7 ]胡小菲, 黄振宇, 洪方文. 螺旋桨非定常力的粘性数值分析[J]. 水动力学研究与进展, 2009, 24(6): 734–739.HU Xiao-fei, HUANG Zhen-yu, HONG Fang-wen. Unsteady hydrodynamics forces of propeller predicted with viscous CFD[J]. Chinese Journal of hydrodynamics, 2009, 24(6):734–739.

[ 8 ]黄振宇. 螺旋桨的非定常力计算[C]//第九届全国水动力学学术会议暨第二十二届全国水动力学研讨会文集. 成都, 2009.HUANG Zhen-yu. Numerical simulation of unsteady force of marine propeller[C]// Proceedings of 9th Hydrodynamics Meeting, 2009.

[ 9 ]刘登成, 洪方文, 张志荣. 等. 伴流中螺旋桨非定常力粘性数值方法研究[C]//第二十三届全国水动力学研讨会暨第十届全国水动力学学术会议文集. 西安, 2011.LIU Deng-cheng, HONG Fang-wen, ZHANG Zhi-rong, et al.Research on viscous numerical method of propeller unsteady force in wake[C]// Proceedings of 10th Hydrodynamics Meeting, 2011.

[10]马超, 肖锋, 张志谊. 等. 艇体非均匀伴流场下螺旋桨压力场[J]. 噪声与振动控制, 2013, 35(2): 49–53.MA Chao, XIAO Feng, ZHANG Zhi-yi, et al. Numerical study on pressure pulsation near the propeller with non-uniform wake field[J]. Journal of Noise and Vibration Control, 2013,35(2): 49–53.

[11]JESSUP S. Measurement of multiple blade rate unsteady propeller force [R]. David Taylor Research Center, Ship Hydromechanics Department Research and Development Report, DTRC-90/015 May 1990.

[12]HUANG T, LIU H L, GROVES N, et al. Measurements of flows over an axisymmetric body with various appendages in a wind tunnel: the DARPA SUBOFF experimental program[C]//Proceeding of 19th Symposium on Naval Hydrodynamics. Seoul, Korea, 1992.

[13]刘志华, 熊鹰, 韩宝玉. 消涡整流片对潜艇马蹄涡的控制及其与辅翼效果的比较[J]. 船舶力学, 2011, 15(10):1102–1105.LIU Zhi-hua, XIONG Ying, HAN Bao-yu. Comparison on the submarine horseshoe vortex control effects by vortex control bafflers and fillets[J]. Journal of Ship Mechanics, 2011, 15(10):1102–1105.