基于CFD的潜艇阻力及流场数值计算

涂海文，孙江龙

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

基于CFD的潜艇阻力及流场数值计算

涂海文，孙江龙

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

运用雷诺平均N－S方程，使用CFD前处理软件ICEM CFD划分流场网格，采用RNG k－ε湍流模型，实现了对裸潜体、带指挥台围壳艇体、带十字尾翼艇体、全附体潜艇4种模型的阻力及粘性流场的数值模拟。通过数值模拟，得到了潜艇表面压力分布情况和附体附近流场的一些特性，为进一步优化潜艇的艇型和分析潜艇的流噪声打下了基础。而阻力的对比在一定程度上验证了数值模拟的可靠性。

潜艇;CFD;RNG k－ε模型;摩擦阻力

0 引言

潜艇周围的流场特别是尾流场的流动特性，不仅对潜艇的水动力性能产生直接的影响，而且引发的流动是潜艇水动力噪声的主要来源之一，对潜艇隐蔽性有重要影响。当航速和排水量一定时，流场的特性与潜艇的主尺度和艇体形状密切相关。因此，开展潜艇几何形状与流场特性的关联研究，寻求能在已知艇体形状特征下精确预报流场特性的数值计算方法，为潜艇概念设计阶段得出阻力小、噪声低的最佳潜艇设计方案提供评价依据，一直是船舶水动力学领域中一个重要的前沿研究课题。

长期以来，经典流体理论无法在工程中直接应用，船舶流体力学工作者主要依靠水池实验来粗略地了解船舶流体性能。但船模实验耗资巨大，实验周期长，且水池实验受到尺度效应、外界干扰等影响，船模实验的弊端显而易见。随着计算机技术的迅猛发展，数值计算方法逐渐成为研究船舶和潜艇流场的另一有效方法［1］。模型试验和CFD数值仿真模拟各有自己特定的优势。现在，科学工作者们往往将CFD仿真结果与模型试验结果相互印证，使这2个方面互相促进和发展。

近年来，潜艇流场的数值模拟得以蓬勃发展。赵峰，周连第［2］运用复杂流场多块耦合计算方法，对由于指挥台而造成的潜艇后体不均匀流场作了较成功的数值模拟。张楠等［4］用数值方法预报了研究用美国DARPA潜艇模型SUBOFF与CSSRC潜艇模型SM－x的阻力与尾流场，并与试验结果进行对比。李艳和姚震球［6］运用数值方法对带附体潜艇尾流场进行模拟与验证。操盛文和吴方良［7］采用不同数量的网格分别对不同尺度的SUBOFF模型在高雷诺数条件下的流场进行数值模拟，研究网格数量和艇体主尺度大小对高雷诺数条件下潜艇阻力计算结果的影响。

本文建立了裸艇体、带指挥台围壳艇体、带十字尾翼艇体、全附体潜艇4个不同的潜艇模型，然后进行CFD的仿真和对比计算。

1 数学模型

1.1 控制方程

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

1.2 湍流模型

上面的方程组是不封闭的，需要采用相应的湍流模型。本文采用RNG κ－ε湍流模型，该湍流模型数学表达可以写成下面的形式［1］:

湍流动能k方程为

在RNG κ－ε模型中，通过在大尺度运动和修正后的粘度项体现小尺度的影响，而使这些小尺度运动有系统地从控制方程中去除。所得到的k方程和ε方程与标准κ－ε模型非常相似。

与标准κ－ε模型比较发现，RNG κ－ε模型修正了湍动粘度，考虑了平均流动中的旋转及旋流流动情况，从而可以更好地处理高应变率及流线弯曲程度较大的流动。

需要注意的是，RNG κ－ε模型仍是针对充分发展的湍流有效，即是高Re数的湍流计算模型，而对近壁区内的流动及Re数较低的流动，必须使用壁面函数法或低Re数的κ－ε模型来模拟。

1.3 数值计算方法

采用有限体积法(FVM)对控制方程进行离散，用SIMPLEC方法计算压力速度关联方程。时间项采用中心差分格式，对流扩散项采用二阶迎风格式。

2 计算模型

2.1 计算对象

建立如图1所示的带全附体的潜艇计算模型［3］。潜艇的指挥台机翼采用的是3:2椭圆首部和NACA0020机翼尾部的复合机翼，则相应的机翼厚度为330 mm，潜艇的十字尾翼采用了NACA0020尾翼。潜艇总长度为995 mm，最大半径为45.6 mm，指挥台机翼距潜艇首部长度为220 mm，机翼后面的平行中体长度为240 mm，指挥台机翼的弦长为140 mm，指挥台高60 mm，十字尾翼最高点距轴中心线55 mm，弦长分布在轴向910～960 mm之间。取潜艇中心线与头部的交点为原点，沿潜艇长度方向定义为x轴正向，机翼附体的高度方向为y轴正向，建立笛卡尔直角坐标系［1］。

图1 全附体潜艇计算模型Fig.1 The computing model of submarine with all appendages

计算域尺寸选取参照图2所示。为了使数值计算结果更接近无限空间内的三维绕流，计算域应取得足够大。为此，上游边界距原点0.3倍艇长，下游边界距原点1.7倍艇长，上、下侧边界距轴线400 mm。

图2 计算域xy平面图Fig.2 Calculation domain of xy plane

光艇体模型、带指挥台围壳艇体模型、带十字尾翼艇体模型均可由全附体模型得到，它们的计算域与全附体计算域相同。

2.2 计算网格

划分网格时主要遵循的原则是:靠近潜艇指挥台、十字尾翼和潜艇表面及交接部等区域采取比较细密的网格，由于壁面及附体附近流场比较复杂，含有分离、旋涡及其脱落等现象，因此在这一区域需要布置较为密集的网格节点，以有利于研究潜艇指挥台交接部和十字尾翼处的流场分析。而在远离壁面及附体的区域，则可布置得较稀疏一些，以节省整个计算的工作量。

网格工具使用的是ICEMCFD网格划分软件。采用非结构性的四面体网格和三棱柱网格，四面体网格适合对结构复杂的几何模型进行快速高效的网格划分，三棱柱网格可以更好地模拟边界层效应。在靠近潜艇指挥台和十字尾翼的区域必须通过网格加密形成更为细密的网格，以便于更好地对这些区域的流场进行模拟。网格图如图3～图5所示。

2.3 边界条件

对于边界条件的定义，为了更好地模拟无限流场的环境，在整个流域中，除了出口，其余2个边界均设置为速度入口。数字计算域的前端面I和侧表面M为来流入口，后端面O为来流出口，见图6。

1)层流边界条件设置

潜艇、指挥台和十字尾翼的壁面设置为固壁条件，满足无滑移条件。

采取层流计算模型时，在非定常计算过程中，每个时步取0.01 s，计算时长为10 s。

2)湍流边界条件设置

计算域来流边界设置为速度入口velocity-inlet，u=U0，v=w=0。其中初始来流速度U0=6.6 m/s，同样以潜艇长作为特征长度，则对应的雷诺数Re=6.6×106。计算域出口处边界条件设置为自由出流outflow。

潜艇、指挥台和十字尾翼设置为壁面边界条件wall，满足无滑移条件，近壁区修正采用加强型壁面函数 Enhanced Wall Treatment。

在非定常计算过程中，每个时步取0.003 s，计算时长为3 s。

3 计算结果分析

3.1 摩擦阻力分析

潜艇运动中所受的总阻力Rt是所有流体作用力沿运动方向的合力亦即潜艇表面上所有微面积ds上切应力τ和压力p在运动方向的合力，即

式中:前1项积分表示由作用在潜艇表面上切应力所造成的阻力，即摩擦阻力Rf;第2项积分表示由作用在潜艇表面上压力所造成的阻力，即粘压阻力Rpv。(τ，x)表示表面剪切应力与 x方向的夹角，(p，x)表示表面压力与x方向的夹角。工程中一般定义摩擦阻力系数及粘压阻力系数如下:

在工程实践中，对于船舶摩擦阻力的预报通常使用“相当平板假定”计算，即利用与船体表面积相当的平板的摩擦阻力来代替船舶的摩擦阻力。

造船界常用的摩擦阻力公式为(ITTC57):

从表1可看出，数值模拟出来的摩擦阻力系数与用经验公式算出来的摩擦阻力系数相比，误差较小，它们的差别在于数值模拟中考虑了潜艇三维船体表面弯曲度的影响，而经验公式未考虑此影响，因此本文的数值模拟结果是合理及可以接受的。从另一方面也验证了数值模拟的可靠性。大，在指挥台围壳和十字尾翼处，压力产生了突变。而对于带指挥台围壳艇体模型，在湍流和层流2种情况下的差别在于湍流的最大压力系数比层流的要大一些，再就是指挥台尾部的“高”压力区域比层流的要小，也就是说湍流情况下压力恢复至潜艇光体相同的压力比层流的要快。从图中还可以看出，4种模型模拟的情况也比较吻合。

1)裸艇体模型静压力系数分布图

?

3.2 壁面静压力系数分析

图7～图14依此描述了裸艇体模型、带指挥台围壳艇体模型、带十字尾翼艇体模型、全附体模型4种模型分别在湍流和层流状态下y=0剖面上下壁面沿x向静压力系数分布情况。从图中可以看出，湍流和层流的计算结果比较相似，均是潜艇头部压力比较

2)带指挥台围壳艇体模型静压力系数分布图

图9 湍流下y=0剖面沿X向静压力系数分布Fig.9 The distribution of static pressure coefficient on y=0 profile along the X direction under the turbulence

图10 层流下y=0剖面沿X向静压力系数分布Fig.10 The distribution of static pressure coefficient on y=0 profile along the X direction under the laminar

3)带十字尾翼艇体静压力系数分布图

4)全附体潜艇模型静压力系数分布图

图13 湍流下y=0剖面沿X向静压力系数分布Fig.13 The distribution of static pressure coefficient on y=0 profile along the X direction under the turbulence

图14 层流下y=0剖面沿X向静压力系数分布Fig.14 The distribution of static pressure coefficient on y=0 profile along the X direction under the laminar

图15和图16将湍流和层流2种情况下全附体潜艇模型静压力系数的分布云图形式清晰，变化的过程很直观。

3.3 流场分析

图17～图19展示了指挥台围壳处的流线分布。从图中可以看出，在围壳头部附近流体出现分流现象，速度方向发生改变，然后在围壳尾部合流，这就是指挥台围壳头部压力很大的原因。

图20展示了十字尾翼在y=0剖面处速度分布云图。图21和图22是十字尾翼在y=0剖面和z=0剖面处的流线图。在y=0剖面u速度分布图上观察到了负速度区，而在流线图21和图22上，表现为明显的涡。

4 结语

1)运用CFD数值模拟技术很好地预报了潜艇的摩擦阻力系数，而通过对潜艇壁面压力的数值预报，可以发现在潜艇的头部，指挥台围壳的头部以及十字尾翼的前端面压力均比较大，从而掌握了潜艇上一些显著地方的壁面压力分布情况，这些都为潜艇快速性和潜艇艇型优化提供了非常好的参考。

2)从层流和湍流2种状态下潜艇的压力分布对比当中可以看出，在2种不同状态下潜艇的压力分布曲线的走势具有一致性，只是压力的大小和一些局部地方的压力分布有所不同。

3)从流场分析中可以看出，在附体的附近流线比较密集，在潜艇的尾翼处有明显的涡现象产生，这就为后面潜艇流噪声的分析作了很好的铺垫。

［1］陈小邹.翼型附体与平板交接部流场分析［J］.舰船科学技术，2010，32(10):8 －11.

CHEN Xiao-zou.Junction flows analysis around flat plate mounted a foilappendage［J］.Ship Science and Technology，2010，32(10):8 －11.

［2］赵峰，周连第.潜艇含指挥台附体区域周围粘性流场的多块耦合计算［J］.水动力学研究与进展(A辑)，1996，11(4):448－458.

ZHAO Feng，ZHOU Lian-di.Multiblock coupled computation of the viscous flow field around a submarine sail zone［J］.Journal of Hydrodynamics，1996，11(4):448 －458.

［3］LEE Y T，HUANG T T，DAVID S C H.Flow prediction about an appended body of revolution［J］.Taylor Research Center Bethesda，MD 20084 －5000.

［4］张楠，沈泓萃，姚惠之.潜艇阻力与流场的数值模型与验证及艇型的数值优化研究［J］.船舶力学，2005，9(1):1－13.

ZHANG Nan，SHEN Hong-cui，YAO Hui-zhi.Validation of numerical simulation on resistance and flow field of submarine and numerical optimization of submarine hull form［J］.Journal of Ship Mechanics，2005，9(1):1 －13.

［5］张楠，沈泓萃，姚惠之，等.带流水孔潜体流场数值模拟［J］.船舶力学，2004，8(1):1 －11.

ZHANG Nan，SHEN Hong-cui，YAO Hui-zhi，et al.Numerical simulation of flow around a submerged body with flood hole［J］.Journal of Ship Mechanics，2004，8(1):1 －11.

［6］李艳，姚震球.带附体潜艇尾流场的数值模拟与验证［J］.江苏科技大学学报(自然科学版)，2006，20(2):7－12.

LI Yan，YAO Zhen-qiu.Numerical simulation and validation of wake field of submarine with stern appendage［J］.Journal of Jiangsu University of Science and Technology(Natural Science Edition)，2006，20(2):7 －12.

［7］操盛文，吴方良.尺度效应对全附体潜艇阻力数值计算结果的影响［J］.中国舰船研究，2009，4(1):33 －37.

CAO Sheng-wen，WU Fang-liang.Investigation of scaling effects on numerical computation of submarine resistance［J］.Chinese Journal of Ship Research，2009，4(1):33 －37.

［8］YANG Chi，LOHNER R.Prediction of flows over an axisymmetric body with appendages［A］.The 8 th International Conferenceon NumericalShip Hydrodynamics，Busan，Korea，September，2003.

［9］李新汶，陈源，王文琦.CFD在潜艇艇附体与艇体连接形式研究中的应用［J］.船舶力学，2003，7(5):28－32.

LI Xin-wen，CHEN Yuan，WANG Wen-qi.Applicationof CFD to junction form of strern appendages.［J］.Journal of Ship Mechanics，2003，7(5):28 －32.

［10］王志博，姚惠之，张楠.指挥台围壳对潜艇尾流影响的计算研究［J］.船舶力学，2009，13(2):196 －202.

WANG Zhi-bo，YAO Hui-zhi，ZHANG Nan.Calculation and analysis of the effects of sail form on submarine wake［J］.Journal of Ship Mechanics，2009，13(2):196 －202.

［11］卢云涛，张怀新，潘徐杰.全附体潜艇的流场和流噪声的数值模拟［J］.振动与冲击，2008，27(9):142 －189.

LU Yun-tao，ZHANG Huai-xin，PAN Xu-jie.Numerical simulation of flow-field and flow-noise of a fully appendage submarine［J］.Journal of Vibration and Shock，2008，27(9):142－189.

Numerical analysis of resistance and flow field of submarine based on CFD

TU Hai-wen，SUN Jiang-long
(College of Naval Architecture and Ocean Engineering，Huazhong University of Science and Technology，Wuhan 430074，China)

The article numerically simulated the resistance and flow field of bare submarine body with different appendages and submarine body full-appended with RNG k － ε model.And solving solutions with RANS equations.The girds were divided by ICEM CFD of CFD pretreatment software.Through numerical simulation，it obtained the submarine surface pressure distribution and some characteristics of the flow field near appendages.This is further to lay the foundation for optimizing the boat type and analyzing noise of submarine that induced by flow.And resistance contrast verified the reliability of numerical simulation to some degree.

submarine;RNG k－ε model;ICEM CFD;frictional resistance

U674.76;TB53

A

1672－7649(2012)03－0019－07

10.3404/j.issn.1672－7649.2012.03.004

2011－05－10;

2011－06－28

涂海文(1986－)，男，硕士研究生，主要从事船舶与海洋工程研究。