载人潜器阻力的数值计算方法分析

李 佳，黄德波，邓 锐

（哈尔滨工程大学多体船技术国防重点学科实验室，哈尔滨 150001）

1 引 言

潜器的种类繁多，用途非常广泛，即使作为同一用途的潜器其三维几何形状也各不相同，而且常常带有附体及支架等复杂结构，因此对潜器水动力性能的研究有一定的困难，针对不同问题使用的研究方法也很多。

计算流体力学（CFD）方法周期短、费用低，是研究潜器水动力性能的有效手段，将CFD方法的计算结果与试验数据进行比较，可以验证CFD方法的有效并实现物理试验所不容易完成的模拟。Le等[1]用数值方法求解了有制导推进器的AUV的定常流动，用有限差分法求解了流函数的轴对称微分方程，并将结果用于计算速度场和压力场。Fuglestad等[2]用CFD方法计算了HUGIN 3000原型的阻力，并与模型试验数据进行了比较，验证了CFD计算的有效性和高雷诺数阻力预报的可能性。张怀新等[3]用有限体积法求解了碟状水下航行器三维粘性绕流场的不可压缩RANS方程，得出了航行器阻力性能的数值计算结果并与试验数据进行了比较。Nishi等[4]用商业CFD软件STAR－CD计算了水下航行器的阻力，并与拖曳水池试验结果进行了比较，考虑阻力和推进的平衡，得出航行器的推进性能。Phillips等[5]用CFD方法计算了三个不同形状和尺寸的AUV的阻力，并将结果与试验数据进行了比较，验证了使用RANS流动求解器完成概念设计的可能性。刘玉秋等[6]采用混合均质流理论、Realizable k-ε湍流模型对某水下非流线型航行体模型进行了自然状态和通气产生超空泡状态下阻力的数值模拟，得到了自然和通气状态下非流线型水下航行体阻力的变化规律。Listak等[7]用CFD方法得出了水下仿生机器人的水动力特性，并将数值模拟得出的阻力值与实际物理测量的阻力值进行了比较。Seo等[8]为了测量水下滑翔物的性能，用CFD处理器得到水动力系数从而对俯仰进行了控制仿真。

本文所研究的潜器是载人潜器，载人潜器在设计上形式多样，关于水动力性能方面的资料在国内外并不多见。本文用数值模拟方法对一艘载人潜器的缩尺模型的水动力性能进行了研究，得到了可行的计算方法和有参考价值的结论。作者首先对潜器计算模型的绕流场进行区域离散，得到适于计算的多块结构化网格。然后采用不同的湍流模型，对潜器周围的绕流场进行数值模拟，并将结果与试验值进行对比。改变入口的湍流强度，探讨入口湍流强度对数值计算结果的影响。

2 区域离散

本文所研究的载人潜器模型的主要尺度参数如表1所示，表中给出了潜器计算模型主体和底部主要附体的主尺度参数，缩尺比为1:6。

三维复杂体的控制域进行区域离散时，生成非结构化网格所用的时间少，而且操作也比较容易，但是往往生成的非结构化网格数量很大，进行数值求解时是一般的PC机或工作站所无法实现的。而且对于生成的非结构化网格很难进行人为控制，所生成的网格会产生局部区域过于密集或系数的现象，网格数量无法预测。

本文所研究的三维载人潜器外型复杂，并且在主体上带有许多附体，因此将控制域离散为非结构网格的数量将会很大。结构网格的数据按照顺序存储，根据数组的下标可以方便地索引和查找。结构网格可以减少相应的存储开销，而且由于网格具有贴体性，流场的计算精度可以大大提高，因此作者花费了大量的时间将控制域离散为结构网格。

潜器周围流场的控制域离散如图1所示，入口边界到潜器最前端距离为L（L为潜器计算模型主体的总长），出口边界到潜器最后端距离为2L，控制域半径为L。由于潜器周围绕流场的数值模拟为不可压缩外流问题，入口边界设置为速度入口，出口边界设置为自由出流，来流方向与控制域轴向平行，周向边界设置为对称边界。

表1 潜器主要尺度参数Tab.1 Main parameters of a manned submersible

控制域采用多块搭接结构化网格进行离散，将控制域根据需要分为96个子块，在不同子块内分别生成结构网格。多块搭接结构网格中相邻两块公共边界上的网格点为两块公共的网格点，不同子块的网格共同充满整个流场，相互之间没有重叠区，相邻块之间以公共边界相连接。在流场计算中，不同块之间的物理信息通过边界点进行交换。作者尽可能描述潜器型体的具体构形，避免过多简化，只在潜器实际几何体的基础上略去较简单的支杆以及较小的推进器。控制域的结构化网格采用O型和H型拓扑结构，潜器主体以及主要附体的表面网格分别见图2～3。

3 基本控制方程

在研究载人潜器的水动力性能时，流体可压缩性的影响是被忽略的。假设流体不可压，则潜器粘性绕流场中流体运动的连续性方程为[9]：

对于不可压缩流体，雷诺平均RANS方程的张量形式为[9]：

其中：ui为速度分量时均值，ui′为速度分量脉动值，为雷诺应力项。

压力和速度用SIMPLEC算法耦合，压力方程离散采用二阶格式。

4 湍流模型

为了得到适用于所研究潜器的数值计算的湍流模型，作者采用了RNG k-ε湍流模型、Realizable k-ε湍流模型、标准k-ω湍流模型以及SST k-ω湍流模型分别对潜器周围绕流场进行数值计算，并根据计算结果改变入口边界的湍流强度从而得到适用于潜器数值计算的方法。

RNG k-ε湍流模型输运方程[10]：

Yakhot等[11]提出重整化群（Renormalization Group）的理论，RNG k-ε模型有效改善了精度，考虑到了湍流漩涡，为湍流Prandtl数提供了一个解析公式。

Realizable k-ε湍流模型输运方程[10]：

Realizable k-ε模型对于平板和圆柱射流的发散比率预测更精确，对于旋转流动、强逆压梯度的边界层流动、流动分离和二次流有很好的表现。

标准k-ω湍流模型输运方程[10]：

Kolmogorov首先提出了k-ω湍流模型的概念，Launder和Spalding对k-ω模型进行了改进，Wilcox进一步对k-ω湍流模型进行应用研究。标准k-ω模型基于Wilcox k-ω模型，是为考虑低雷诺数、可压缩性和剪切流传播而修改的。标准k-ω湍流模型能够很好地处理近壁面流动问题，但对自由来流有很大的依赖性。

SST k-ω湍流模型输运方程[10]：

SST k-ω模型是由Menter发展的，考虑了湍流剪应力的输运特性，能够准确地预报由于逆压梯度导致的流动分离点和分离区域。

5 自由来流湍流度

如果两种流体运动具有几何相似的边界，并且它们的速度场是几何相似的，也就是说，如果它们有几何相似的流线，则这两种流体运动是动力学相似的。根据动力学相似原理，如果两个物体几何相似，其它相关条件如相对粗糙度、自由来流湍流度、空泡数等相同，则物体的阻力系数就可表示为：

上式中的四个变量中可以略去次要因素而考虑主要因素来实现流动相似，对于载人潜器阻力性能的研究，我们考虑的主要因素是雷诺数Re。

当动力学相似所考虑的其它条件相同时，物体在静水中（湍流度为0）航行和自由来流流经静止物体时物体的阻力系数是不同的[12]，这是因为自由来流的湍流度不同。也就是说载人潜器的计算模型的数值模拟无法做到与实际情况相同，自由来流的湍流度不能设置为0，因此入口边界的湍流度的取值会对模拟结果产生影响。

充分发展的管内流动可以根据半经验公式给出入口边界处的湍流强度和湍流长度尺度，而不可压缩外部绕流只能根据经验或试验给出入口的湍流强度I和湍动粘度比μt/μ，湍流强度应为很小的量，而湍动粘度比应在1～10之间。为了研究适于潜器流场数值计算的入口边界条件，入口处湍流强度和湍动粘度比的设置首先根据文献[13]取I=0.5%，μt/μ=5，然后改变入口处的湍流强度I，用不同的湍流模型进行潜器周围流场的数值计算，将计算结果与试验值相比较，得出适于计算的入口湍流强度。

6 计算结果及分析

为了探讨自由来流湍流强度对潜器阻力数值模拟结果的影响，以及不同湍流模型对自由来流湍流度的依赖程度，首先设置入口湍流强度I=0.5%以及I=0.1%，分别采用RNG k-ε湍流模型、Realizable k-ε湍流模型、标准k-ω湍流模型以及SST k-ω湍流模型对潜器周围绕流场进行数值模拟。阻力计算结果如图4～5所示，其中Vm为来流速度，Rt为潜器总阻力。对比可见，改变入口湍流强度时，采用标准k-ω湍流模型得到的潜器阻力计算结果差别很大，而采用其它三种湍流模型得到的计算结果基本没有差别。下面对标准k-ω湍流模型对自由来流湍流强度的敏感程度作简要分析。

根据上面的计算结果，入口湍流强度的变化范围取为I=0.1%-0.5%，不同来流速度下采用适当的自由来流湍流强度，用标准k-ω湍流模型得出的潜器阻力计算结果如图6所示。相应速度下设置的入口湍流强度与雷诺数的关系曲线见图7，其中Re=VmL/ν为雷诺数，运动粘性系数ν=1.12434×10-6m2/s，图中还给出了该关系曲线的拟合趋势曲线。由图6可知，设置适当的入口湍流强度采用标准k-ω湍流模型对潜器模型的阻力计算的结果与试验值是非常吻合的。对于不可压缩外流而言，入口湍流强度的取值并不像内流一样有经验公式可循，我们只能根据试验值找出一定的规律。图7中自由来流湍流强度取值的趋势曲线的函数为I=2×107×Re-1.2298（%），此函数关系是对比载人潜器模型的阻力试验值，采用标准k-ω湍流模型进行载人潜器模型的阻力计算时对于不同来流速度下入口湍流强度的设置依据。

采用RNG k-ε湍流模型、Realizable k-ε湍流模型以及SST k-ω湍流模型对潜器的阻力进行数值计算时，自由来流的湍流强度分别取I=0.5%和I=0.1%，阻力的计算结果变化很小。因此自由来流湍流强度的变化范围在I=0.1%-0.5%时，采用这三种湍流模型对潜器进行阻力计算的结果将不会有过大差别。

图8～9给出了潜器模型速度为Vm＝2.6m/s，采用Realizable k-ε湍流模型，入口湍流强度设置为I＝0.5%时潜器主体和主要附体的压力云图，由图中可以看到附体前后压差很大，并且附体之间的流场产生了相互干扰。

7 结 论

为了得到适于潜器周围绕流场进行数值模拟的计算方法，本文对某载人潜器周围控制域进行多块搭接结构化网格离散，分别采用RNG k-ε湍流模型、Realizable k-ε湍流模型、标准k-ω湍流模型以及SST k-ω湍流模型对潜器周围绕流场进行数值模拟，并改变自由来流的湍流强度的设置，比较了阻力计算结果并将其与试验值进行了对比。计算结果分析表明：

（1）采用标准k-ω湍流模型对潜器周围绕流场进行数值计算时，自由来流湍流强度在I=0.1%-0.5%范围内变化，得到的潜器阻力计算结果有所不同。根据阻力试验值，可以给出自由来流湍流强度I=2×107×Re-1.2298（%），采用标准k-ω湍流模型进行潜器阻力的数值计算时不同速度下入口湍流强度的设置要有变化，可以参照此函数关系进行设置。

（2）采用RNG k-ε湍流模型、Realizable k-ε湍流模型以及SST k-ω湍流模型对潜器的阻力进行数值计算时，自由来流湍流强度在I=0.1%-0.5%范围内变化对潜器阻力的计算结果影响很小。采用这三种湍流模型进行潜器阻力的数值计算时，可以将入口湍流强度设置为I=0.5%。

本文对采用标准k-ω湍流模型进行潜器阻力数值计算情况时，是在有试验值可依的情况下得出的自由来流湍流强度的设置随雷诺数的变化规律，此规律还要经过大量的计算以及试验验证。

[1]Le Page Y G,Holappa K W.Hydrodynamics of an autonomous underwater vehicle equipped with a vectored thruster[C]//Oceans Conference Record(IEEE).Oceans 2000,Providence RI USA,2000:2135-2140.

[2]Fuglestad A L,Grahl-Madsen M.Computational fluid dynamics applied on an autonomous underwater vehicle[C]//Proceedings of the 23rd International Conference on Offshore Mechanics and Arctic Engineering.Vancouver BC Canada,2004,USA:American Society of Mechanical Engineers,2004:447-451.

[3]张怀新,潘雨村.圆碟形潜水器阻力性能研究[J].上海交通大学学报,2006,40(6):978-982,987.Zhang Huaixin,Pan Yucun.The resistance performance of a dish-shaped underwater vehicle[J].Journal of Shanghai Jiaotong University,2006,40(6):978-982,987.

[4]Nishi Y,Kashiwagi M,Koterayama W,et al.Resistance and propulsion performance of an underwater vehicle estimated by a CFD method and experiment[C]//Proceedings of the 17th 2007 International Offshore and Polar Engineering Conference.Lisbon Portugal,2007,Cupertino,CA 95015-0189,United States:International Society of Offshore and Polar Engineers,2007:2045-2052.

[5]Phillips A,Furlong M,Turnock S R.The use of computational fluid dynamics to assess the hull resistance of concept autonomous underwater vehicles[C]//OCEANS 2007-Europe.Aberdeen Scotland United kingdom,2007.

[6]刘玉秋,于开平,张嘉钟.水下非流线型航行体减阻的数值模拟与比较[J].工程力学,2007,24(2):178-182.Liu Yuqiu,Yu Kaiping,Zhang Jiazhong.Numerical simulation and comparison on drag reduction of non-streamline submerged body[J].Engineering Mechanics,2007,24(2):178-182.

[7]Listak M,Pugal D,Kruusmaa M.CFD simulations and real world measurements of drag of biologically inspired underwater robot[C]//US/EU-Baltic International Symposium:Ocean Observations,Ecosystem-Based Management and Forecasting.Tallinn Estonia,2008.

[8]Seo D C,Jo G,Choi H S.Pitching control simulations of an underwater glider using CFD analysis[C]//Oceans’08 MTS/IEEE Kobe Technology Ocean’08 Voyage toward the Future.Kobe Japan,2008.

[9]王福军.计算流体动力学分析—CFD软件原理与应用[M].北京:清华大学出版社,2004.Wang Fujun.Computational fluid dynamics analysis－principle and application of CFD software[M].Beijing:Tsinghua U-niversity Press,2004.

[10]Fluent Inc.,FLUENT User’s Guide[K].Fluent Inc,2003.

[11]Yakhot V,Orzag S A.Renormalization group analysis of turbulence.Ⅰ.Basic theory[J].Journal of Science Computing,1986,1(1):3-51.

[12]Douglas J F,Gasiorek J M,Swaffield J A.Fluid Mechanics[M].Great Britain:The Pitman Press,1979.

[13]韩占忠,王 敬,兰小平.FLUENT—流体工程仿真技术与实例应用[M].北京:北京理工大学出版社,2004.Han Zhanzhong,Wang Jing,Lan Xiaoping.FLUENT－Emulational technique and instance application of fluid engnieering[M].Beijing:Beijing Insitute of Technology Press,2004.