基于CFD的潜艇模型自航仿真分析

杨 琴，王国栋，张志国，冯大奎，王先洲

华中科技大学船舶与海洋工程学院，湖北武汉430074

0 引 言

对潜艇周围的流场信息特别是尾部流动特性进行研究一直是潜艇设计和流体力学领域的一项重要课题。当潜艇尾部螺旋桨盘面处的伴流场不均匀和不稳定时，会降低螺旋桨的推进效率。研究潜艇螺旋桨与尾部操纵面之间的互相影响，对提高潜艇的推进效率和机动性具有重要意义。可通过潜艇的自航试验来获得螺旋桨与船体之间的相互作用等诸因素，这是分析研究潜艇快速性与推进效率的重要手段。

数值水池的重要基础是CFD 技术，即应用CFD 方法进行船舶流动数值模拟来实现船舶水动力学性能的数值计算与预报。数值水池船舶阻力与螺旋桨敞水性能计算发展较早。对全附体潜艇的数值模拟一般是采用RANS 方程结合k-ε 湍流模型求解［1-2］，或采用大涡模拟方法。目前，螺旋桨的敞水性能研究已发展得较为成熟，一般采用混合网格［3］来对大侧斜螺旋桨进行数值计算，常用的湍流模型为k-ε，SST 和雷诺应力模型［4-5］，大涡湍流模型［6］也被运用到了大侧斜螺旋桨的数值计算中，以上数值计算结果与实验室值都已达到了较高的精度［7-8］。数值水池自航试验研究需实现螺旋桨与船体的整体求解，Choi 等［9］分别采用体力桨和滑移网格方法对船—桨—舵之间的水动力相互影响进行了计算，并与模型试验值进行了对比，获得了比较准确的计算结果；刘祥珺等［10］通过求解RANS 方程，利用VOF 方法追踪自由液面，在数值水池实现了船模的自航试验研究。螺旋桨旋转域与静止域采用滑移网格技术。在既定航速下，船体阻力和螺旋桨推力为螺旋桨转速的函数，通过变化转速，可以得到自航点。采用数值方法得到的自航转速与试验值吻合良好，张楠等［11］采用数值模拟方法研究了潜艇近水面航行时的艇/桨干扰特性。流场采用RANS 方法结合RNG k-ε 湍流模型求解；自由液面捕捉采用VOF方法；螺旋桨采用滑移网格方法。

1 理论基础

本文采用RNG k-ε 湍流模型求解Reynolds平均Navier-Stokes 方程，建立了相关的仿真数学模型，对带七叶螺旋桨的SUBOFF 全附体模型的三维粘性流场进行了数值分析，并采用准静态的方法对潜艇模型的阻力和螺旋桨推力进行了耦合，得到了指定速度所对应的螺旋桨推力和转速。

1.1 控制方程

本文研究分析的“艇—桨”自航条件下螺旋桨的控制方程为不可压缩牛顿流体流动的连续性方程和RANS 方程：

1.2 湍流模型

在涡粘模型中，根据Boussinesq 提出的涡粘假定，引入了湍动粘度，该假定建立了雷诺应力与平均速度梯度的关系，即

通常，涡粘模型有零方程模型、一方程模型和二方程模型。本文采用的RNG k-ε 湍流模型对应的方程组如下：

1.3 敞水特征曲线

螺旋桨模型单独在均匀水流中试验称为敞水试验，由螺旋桨敞水试验得到的是螺旋桨的推力系数KT、扭矩系数KQ和敞水效率η0相对于进速系数J 的变化规律，即螺旋桨敞水特征曲线。与敞水试验对应的船舶计算流体力学（CFD）计算工作称之为螺旋桨敞水性能计算。采用CFD 方法计算螺旋桨的水动力性能，所采用的计算模型具有较高的可靠性，其结果也具有一定的准确性。

螺旋桨的敞水性能包括螺旋桨的推力系数KT、扭矩系数KQ和敞水效率η0，具体计算公式如下：

式中，

式中：J 为进速系数，是影响螺旋桨性能的重要参数，其相当于机翼理论中攻角的概念；Va代表螺旋桨的进速；T 为螺旋桨产生的轴向推力；Q 为螺旋桨运行中产生的扭矩；ρ 为工作介质（20 ℃水）的密度；D 为螺旋桨直径；n 为螺旋桨转速。

1.4 螺旋桨与艇的相互作用

伴流的大小通常用伴流速度u 与船速V 的比值ω 来表示。ω 称为伴流分数，即

螺旋桨在船后工作时引起的船体附加阻力称为阻力增额，螺旋桨发出的推力中只有一部分是用于克服阻力R 并推船前进的，称这一部分力为有效推力T 。通常，将阻力增额称为推力减额，并以ΔT 表示。推力减额ΔT 与推力T 的比值称为推力减额分数：

2 计算对象

本文的计算对象为带桨SUBOFF 全附体模型。该模型艇长4.356 1 m，螺旋桨采用改进的IN⁃SEANE 1619 型七叶螺旋桨。为配合全附体潜艇的尾端直径，建立了直径为302.6 mm 的七叶螺旋桨三维模型，其主要几何参数如表1 所示。

表1 螺旋桨几何参数Tab.1 Geometric parameters of propeller

本文采用SolidWorks 软件对计算对象进行了三维建模，所建立的带桨SUBOFF 全附体三维模型示意图如图1 所示。图中，x 轴与螺旋桨旋转轴一致，以指向船尾为正。

图1 带桨SUBOFF 三维模型Fig.1 3D model of SUBOFF with propeller

3 数值模拟方法

3.1 计算域与网格划分

为了模拟带桨艇周围的流场，进行数值计算时要确定一个大小合适的流场控制体。本文将带桨艇模放在了一个与螺旋桨同轴线的圆柱形流场区域中，进流面取为上游3 m，出流面取为下游5 m，径向取为2.5 m。其中，螺旋桨被一圆柱小体包裹，此为旋转区域，长0.18 m，直径0.5 m，如图2所示。

计算采用结构化网格与非结构化网格相结合的方式，在包裹螺旋桨的小体内采用非结构化网格，其余部分采用结构化网格。进行了合理的布置并选择了适当的网格尺度，以期取得很好的计算精度和计算时间。图3 给出了带桨潜艇表面的网格分布，螺旋桨附近非结构化网格数均为50 万左右，整个计算域的网格数约为250 万。

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

图3 带桨潜艇表面的网格划分Fig.3 Surface mesh on the SUBOFF simulation model

3.2 边界条件

本文对SUBOFF主艇体的数值模拟边界如下：

入口条件（Velocity Inlet）：入口及计算域周向设置为速度入口，速度大小为Va= 8 kn，速度方向平行于桨轴，为x 轴正方向。

出口条件（Pressure Outlet）：认为流动在该处已充分发展，故边界条件设置为压力出口。

壁面条件（Wall）：螺旋桨表面和艇体外均设置为无滑移壁面条件。

Interface：由于使用滑移网格模拟螺旋桨的旋转运动，故螺旋桨旋转区域与外流场区域交界面均设置为interface。

坐标系：坐标原点在桨盘面圆心处，x 轴沿来流正向布置，坐标系符合右手定则。

3.3 自航点的确定

在已知不同航速下艇体阻力曲线和螺旋桨敞水性能特征的基础上，对“艇—桨”整体进行数值模拟的自航试验。自航点的寻找流程如图4所示。

具体步骤如下：

1）由阻力曲线查得给定航速Va= 8 kn 下的艇体阻力Ds。

2）依据经验给定初始推力减额系数t0=0.079，由全附体阻力Ds和t0求出螺旋桨初始推力T0，然后再由给定螺旋桨敞水性能曲线查得效率最高点对应的推力系数KT，最后，由T0和KT求得初始转速n0。

图4 确定自航点流程示意图Fig.4 Flow chart of finding out self-propulsion point

3）给定Va和n0后，对“艇—桨”整体进行RANS 模拟。待迭代收敛后，读取艇体阻力D 和螺旋桨推力T 值。当T ＞D 时，减小转速n ；反之，则增加转速n。改变转速后，再次以初始计算流场为初值进行RANS 模拟，并重复该过程进行迭代计算，直至推力与阻力的平均值相等或达到允许的不平衡度为止。

4）找到自航点后，求取螺旋桨旋转域进流面处的有效伴流系数ωe、螺旋桨推力减额系数t 和相对旋转效率ηR。

4 数值计算结果与分析

4.1 艇带桨水动力性能预报

对螺旋桨的定常性能模拟采用FLUENT 提供的滑移网格模型。采用RNG k-ε 湍流模型对SUBOFF 全附体模型进行阻力计算。在以上基础上，对带桨潜艇在来流速度Va= 8 kn 下进行数值模型的自航试验，取螺旋桨转速n = 13.2 为初始转速，并预报其水动力性能。

本文以DTMB 4119 等螺旋桨为对象进行了研究。采用本文的数值模拟方法计算所得的结果与已知实验值［12］的对比如表2 所示。SUBOFF 全附体模型阻力的数值模拟结果与试验数据的比较如图5 所示。以上误差均控制在3%以内，充分证明了本文采用的数值模拟方法具有较高的准确性。

图5 全附体Suboff实验值与数值计算结果对比Fig.5 Comparison between experiment and numerical simulation of SUBOFF

在以上基础上，对SUBOFF 全附体螺旋桨模型的计算网格，将其Y+值控制在了3.03～41.92 范围内。图6 给出了螺旋桨压力面与吸力面上的压力分布云图。从图中可以看出，7 个桨叶的压力分布基本相同，且压力面从随边至导边，从叶根至叶梢，压力是逐渐增加的，叶根处的压力最小；在吸力面，压力是中间低四周高，靠近叶梢的部分存在低压区。对于一个叶切面而言，压差最大点位于切面最大厚度处。由图7 给出的带桨潜艇表面压力分布图可看出，在艇前、围壳前和机翼前均有较大的压力梯度。由尾流处的轴向速度分布（图8）可以看到，螺旋桨附近有较大的速度变化。而由于螺旋桨的旋转运动，桨后尾流速度场呈绕轴线螺旋状分布，其尾流流线如图9 所示。

图6 螺旋桨吸力面和压力面压力分布云图Fig.6 Pressure distribution on the blade surface of suction side and pressure side

图7 带桨潜艇表面压力分布图Fig.7 Pressure distribution on the surface of SUBOFF

图8 尾流处轴向速度分布图Fig.8 Axial velocity distribution

图9 带桨艇流线图Fig.9 Velocity streamlines

4.2 数值自航计算与结果分析

由给定的潜艇阻力曲线，得到对应的来流速度Va= 8 kn 下SUBOFF 全附体阻力Ds= 185 N。依据潜艇推力减额分数设计参考值（0.1～0.18）［13］，设定初始推力减额分数t = 0.14。在给定螺旋桨的敞水性能曲线上，查得最大效率点取KT值，得到初始转速n = 0.221 7 r/min。根据自航点确定流程，得到带桨的潜艇阻力与螺旋桨推力之间的关系如图10 所示。

图10 数值自航结果分析Fig.10 The analysis of numerical self-propulsion result

由图10 可知，在推力与阻力相等时，即为潜艇自航点，螺旋桨转速n=0.220 7 r/min(J=0.998)，潜艇阻力与螺旋桨的推力值均为213 N。求得螺旋桨旋转域进流面处的有效伴流系数ωe= 0.17，螺旋桨推力减额系数t = 0.13，相对旋转效率ηR=91%，根据数值模拟得到的上述结果，与常规设计手册提供的参数选择范围进行比较，均在合理的取值范围内。

5 结 语

本文采用数值模拟方法系统地研究了全附体潜艇+螺旋桨的水动力特性，清晰、形象地描述了带桨潜艇表面的压力分布情况，以及桨后尾流速度场的绕轴线螺旋状分布规律。对螺旋桨在艇后非均匀进流条件下的推力和力矩的脉动特性进行了分析，得到了可靠的带桨潜艇阻力和螺旋桨推力数据，潜艇自航点在螺旋桨转速n = 0.220 7 r/min（J = 0.998）时，潜艇阻力与螺旋桨的推力值均为213 N。根据潜艇数值自航试验结果，可以通过积分的方法准确计算伴流分数与推力减额分数，详细分析螺旋桨叶片表面的压力分布，梢涡特性和压力脉动特性。该方法的开发可以显著减少试验研究的工作量，缩短研究周期，具有广阔的应用前景。

［1］吴方良，吴晓光，马运义，等.基于雷诺应力方程模型的全附体潜艇尾部伴流场三维粘性数值计算［J］.中国舰船研究，2007，2（6）：4-8.WU Fangliang，WU Xiaoguang，MA Yunyi，et al. Nu⁃merical calculation of 3D viscous wake field of subma⁃rine with full appendages based on Reynolds stress equation model［J］. Chinese Journal of Ship Research，2007，2（6）：4-8.

［2］吴方良，吴晓光，马运义，等.潜艇实艇阻力预报方法研究［J］.中国舰船研究，2009，4（3）：28-32.WU Fangliang，WU Xiaoguang，MA Yunyi，et al. Meth⁃od of predicting the total submerged resistance of sub⁃marines［J］. Chinese Journal of Journal of Ship Re⁃search，2009，4（3）：28-32.

［3］RHEE S H，LEE C，LEE H B，et al. Rudder gap cavita⁃tion：Fundamental understanding and its suppression devices［J］. International Journal of Heat and Fluid Flow，2010，31（4）：640-650.

［4］杨琼方，王永生，张志宏，等.大侧斜螺旋桨敞水性能的RANSEs 模拟［J］. 湖南科技大学学报（自然科学版），2008，23（4）：66-70.YANG Qiongfang，WANG Yongsheng，ZHANG Zhi⁃hong，et al. RANSEs simulation of highly skewed pro⁃peller’s open water performances［J］. Journal of Hu⁃nan University of Science and Technology：Natural Sci⁃ence Edition，2008，23（4）：66-70.

［5］赵鹏飞. 船用螺旋桨敞水和空化性能CFD 预报［D］.大连：大连理工大学，2011.ZHAO Pengfei. CFD predictions of open water and cav⁃itation performances of marine propellers［D］. Dalian ：Dalian University of Technology，2011.

［6］FELICE F D，FELLI M，LIEFVENDAHL M，et al.Numerical and experimental analysis of the wake be⁃havior of a generic submarine propeller［C］//First Inter⁃national Symposium on Marine Propulsors（SMP'09）.Trondheim，Norway，2009.

［7］ANDERSEN P，KAPPEL J J，SPANGENBERG E.As⁃pects of Propeller developments for a submarine［C］//First International Symposium on Marine Propulsors（SMP'09）.Trondheim，Norway，2009.

［8］胡磊，谭廷寿.大侧斜螺旋桨敞水性能分析［J］.船海工程，2010，39（1）：41-44.HU Lei，TAN Tingshou. Analysis of open water perfor⁃mance of the high-skew propeller［J］. Ship and Ocean Engineering，2010，39（1）：41-44.

［9］CHOI J E，MIN K S，KIM J H，et al. Resistance and propulsion characteristics of various commercial ships based on CFD results［J］. Ocean Engineering，2010，37（7）：549-566.

［10］刘祥珺，孙存楼. 数值水池船模自航试验方法研究［J］.舰船科学技术，2011，33（2）：28-31.LIU Xiangjun，SUN Cunlou. Research on ship self-propulsion model test in numerical tank［J］. Ship Science and Technology，2011，33（2）：28-31.

［11］张楠，张胜利，沈泓萃，等.带自由液面的艇/桨干扰特性数值模拟与验证研究［J］. 水动力学研究与进展A 辑，2012，27（1）：94-99.ZHANG Nan，ZHANG Shengli，SHEN Hongcui，et al.Numerical simulation of hull/propeller interaction with free surface［J］.Chinese Journal of Hydrodynam⁃ics，2012，27（1）：94-99.

［12］盛振邦，刘应中. 船舶原理［M］. 上海：上海交通大学出版社，2004.

［13］BOSWELL R J. Design，cavitation performance，and open water performance of a series of research skewed propellers. NSRDC Report No. 3339 ［R］. 1971：4381-4383.