潜艇阻力及流场数值仿真策略优化分析

图1 几何模型及计算域

Fig.1 Geometry model and computational domain

2.2 网格划分和工况设定

为开展后续的收敛性分析,定义并划分4套网格方案,网格尺寸以 $\sqrt[]{2}$ 为比例逐渐缩小,4套网格方案的具体信息如表1所示,各网格方案的y⁺均小于1^[13],边界层厚度依次为1,1.5,2及 2.5 μm,边界层均布置为21层.图2所示为网格方案3的网格简图,需要提及的是由于计算机软硬件产生的摩尔纹,图2中网格分辨率降低且出现异常纹路, 但数值计算过程中网格分辨率未受影响.除开展相应的网格收敛性计算,后续同样开展时间步长收敛性计算,时间步长同样采用4套方案,分别为 0.0001,0.0005,0.0010及 0.0050 s.

表1 4套网格方案的网格信息

Tab.1 Grid information of four grid schemes

网格方案	围壳表面/个	尾舵表面 (单个舵)/个	艇体表面/个	计算域/ 个
1	9800	2600	267000	6930000
2	11600	3100	294000	8130000
3	13000	4100	311000	9820000
4	15100	5300	352000	12000000

图2

图2 网格方案3简图

Fig.2 Sketch of grid scheme 3

为统一和方便作者进行数值和试验研究,结合计算和实际条件需求,对航速为2.53 m/s的潜艇模型的水动力性能及流场特征进行收敛性计算,以对比不同湍流模型获得的仿真结果.

3 收敛性分析

开展空间和时间敏感性分析是降低数值仿真离散误差最直接的手段,即仿真结果的真值不依赖于空间和时间分辨率,故网格收敛性研究应包含在数值仿真的过程中^[14,15].本节对Sub-off模型航速为2.53 m/s工况下的水动力性能进行收敛性计算,计算使用的是LES-Samgorinsky湍流模型,同时定义总阻力系数 $C_{R} = R / (0.5 ρ U_{ref}^{2} L^{2})$ ,黏性阻力系数 $C_{R V} = R_{V} / (0.5 ρ U_{ref}^{2} L^{2})$ ,压差阻力系数 $C_{Rp} = R_{p} / (0.5 ρ U_{ref}^{2} L^{2})$ ,其中 $R$ 为艇体总阻力, $R_{V}$ 为艇体黏性阻力, $R_{p}$ 为艇体压差阻力,并以上述3个无因次系数作为收敛性分析的评判系数.4套网格的计算结果如表2所示.

表2 4套网格方案的计算结果

Tab.2 Calculation results of four grid schemes

评判系数	方案1/个	方案2/个	方案3/个	方案4/个
1000C_R	1.133	1.156	1.169	1.175
1000C_R	0.962	0.974	0.984	0.989
1000C_R	0.171	0.182	0.185	0.186

网格收敛率 $K = ε_{i, i + 1} / ε_{i, i - 1}, ε_{i, i + 1}$ 是较密两方案结果差, $ε_{i, i - 1}$ 是较稀两方案结果差. $K$ 的计算结果可能是 $K > 1,0 < K < 1$ 或 $K < 0^{}$ ^[14,16-17].网格收敛性计算结果如表3所示.从表3中的计算结果可知选用的网格呈现单调收敛,采用理查森外推法^[16]评估网格的收敛阶P.首先须评估离散阶:

$P = \ln (ε_{i, i + 1} / ε_{i, i - 1}) \ln δ$

式中:δ是网格缩比,取 $\sqrt[]{2}$ .

表3 网格收敛性计算结果

Tab.3 Calculation results of grid convergence

网格方案	参数	1000C_R	1000C_R_V	1000C_Rp
1~3	K	0.565	0.833	0.273
	P	1.646	0.526	3.749
	I	0.021	0.063	0.001
2~4	K	0.462	0.500	0.333
	P	2.231	2.000	3.170
	I	0.006	0.006	0.001

评估3套网格的收敛指数:

$I_{i, i + 1} = D_{S} |e_{i, i + 1}| / (δ^{P} - 1)$

式中:D_S=1.25为网格收敛安全系数^[18]; $e_{i, i + 1} = ε_{i, i + 1}$ .理论上,当 $P = 2$ 时,网格正交性或湍流模型非线性对计算结果准确性的影响程度最小^[14,18],选用中等密度网格进行数值计算具有较高可靠性.可知方案2~4具有更好的收敛率且离散精度更接近理论值,故后续的计算中采用网格方案3进行计算.

表4是不同时间步长计算得到的判定系数结果,以0.001 s的计算结果为基准进行对比.可以看出, 不同时间步长对计算结果的影响较小,为提高计算效率,后续的计算中均采用0.001 s为时间步长.

表4 不同时间步长计算结果

Tab.4 Calculation results of different time steps

时间步长/s	1000C_R	1000C_R_V	1000C_Rp
0.0001	1.162	0.982	0.180
0.0005	1.164	0.983	0.181
0.0010	1.169	0.984	0.185
0.0050	1.175	0.987	0.188

4 结果与讨论

4.1 水动力性能

不同湍流模型对近壁面流动、黏性应力等的解析不同,造成仿真得到的阻力有所差异,表5为不同湍流模型计算得到的潜艇阻力.

表5 不同湍流模型的艇体阻力仿真结果

Tab.5 Simulation results of submarine resistance using different turbulence models

湍流模型	总阻力/N	黏性阻力/N	黏性阻力占总阻力比值/%	压差阻力/N	压差阻力占总阻力比值/%
LES-S	33.617	5.322	15.8	28.295	84.2
LES-W	11.727	4.551	38.8	7.716	61.2
LES-DS	12.311	4.914	39.9	7.397	60.1
DES-SA	34.411	4.783	13.9	29.628	86.1
DES-EBKE	32.632	4.736	14.5	27.896	85.5
DDES	31.133	4.827	15.5	26.306	84.5
IDDES	30.328	4.843	16.0	25.485	84.0
RANS-k-ω	34.403	5.138	14.9	29.265	85.1
RANS-k-ε	31.137	4.772	15.3	26.365	84.7
RANS-SA	35.265	4.899	13.9	30.366	86.1

可见,LES下的两个模型对阻力的仿真结果明显与其他模型不同,两个模型对黏性阻力的估计明显偏大,而对于总阻力的仿真结果更小于其他模型.文献[16]中测得同一弗劳德数 $Fr = U_{ref} / \sqrt[]{gL}$ (g为重力加速度)下的 $C_{R} = 1000 R / (0.5 ρ U_{ref}^{2} L^{2})$ 为1.159.表6为本文不同湍流模型的结果和实验值的对比,表中 $C_{R S}$ 是本文模拟结果.对比发现:LES-W、LES-DS的仿真结果明显偏离实验结果,DDES、IDDES、RANS-k-ε及RANS-SA的计算结果同实验值相比超过5%,高于通常工程应用或学术研究的要求.值得注意的是,LES-S的仿真结果和实验结果最为接近,具备良好的预报精度.

表6 不同湍流模型计算结果与实验结果对比

Tab.6 Comparison of calculation results and experimental results of different turbulence models

湍流模型	C_R	C_R_S	误差/%
LES-S	1.159	1.170	0.9
LES-W	1.159	0.428	63.0
LES-DS	1.159	0.408	64.8
DES-SA	1.159	1.198	3.3
DES-EBKE	1.159	1.136	2.0
DDES	1.159	1.083	6.5
IDDES	1.159	1.055	8.9
RANS-k-ω	1.159	1.197	3.3
RANS-k-ε	1.159	1.084	6.5
RANS-SA	1.159	1.227	5.9

4.2 流场特征

采用Q标准来判定流体为涡旋或应变主导^[19],图3所示为不同湍流模型仿真计算得到的Sub-off周围的涡场.图3清晰地表明LES-W、LES-DS对潜艇重要涡特征——马蹄涡的捕捉不足,不同的是除此外的8个湍流模型均很好地捕捉到了围壳附近马蹄涡的流动特征,另外这两种湍流模型预报的潜艇尾肩部的涡结构亦没有其他模型强烈.选用的10种湍流模型中,RANS类湍流模型对潜艇尾部涡强度的预报结果显然低于其他模型,这可能是RANS速度平均后的结果,是RANS模型的算法特点所决定.值得注意的是,10个模型均没能很好地捕捉或预报围壳稍涡的演化,但对于尾舵稍涡的演化均能较好地捕捉.推测原因:一是围壳稍涡的强度没有尾舵稍涡强度大,二是空间分辨率或时间分辨率不足以捕捉其演化过程.限于计算资源以及本文研究重点,在此不做深入探讨.

图3

图3 不同湍流模型计算艇体周围Q采样体积(Q=5 s^-2)

Fig.3 Calculation of Q-sampling volume around submarine of different turbulence models

伴流是考核潜艇综合性能的重要指标,对伴流的准确预报十分重要.图4所示为不同湍流模型计算得到的Sub-off尾流伴流分数δ随距离变化的曲线,图中,0为潜艇中轴线的最尾端,D、2D、3D、4D为沿着x正方向与潜艇中轴线最尾端的距离.可见,LES-DS、LES-W的预报结果与其他模型有明显偏差,但同属大涡模拟类中LES-S的预报结果有良好的“随群”性,结合对于宏观力的预报结果可知,LES-S模型在预报潜艇宏观变量方面有着较好的适用性.

图4

图4 不同湍流模型预报尾流伴流(0位于艇体尾端)

Fig.4 Wake predicted by different turbulence models (0 at stern of submarine)

为了深入探究不同湍流模型对Sub-off流场的预报结果, 本文提取了如图5虚线所示位置的流场时均变量并加以分析,时均流场的采样总长是流经一倍艇长时间. 图6所示为时均轴向速度的周向变化曲线,图中θ为图5中提取截面沿着白色虚线环向位置的角度,u为流场局部轴向速度,通过u/U_ref对轴向速度进行无因次化对比, 分析其变化规律,为了有效说明流场,将Sub-off周围流场分为内流场(r/D<0.501)和外流场(r/D>0.501),其中r为对应曲线的半径值.可知,10种湍流模型对Sub-off的外流场及其平行中体内流场的预报结果吻合程度非常高,推测由于远离壁面或壁面流动简单,各湍流模型对其处理无明显差异,但即使对平行中体内流场的流动预报差异微弱,Sub-off尾内流场的预报差异却十分明显.这表明,湍流模型对因为几何特征变化引发流动变化更为敏感,对流场“传递过程”引发的变化敏感性较弱.图7所示为流场时均压力的周向分布曲线,可以看出,Sub-off周围流场的时均压力系数 $C_{p} (C_{p} = p / (0.5 ρ U_{ref}))$ 分布曲线对湍流模型地变化不敏感,各个模型均有良好的“随群”性,由于尾部附体以及尾部几何收缩造成的负压力梯度,在x/L=0.9处,流场压力波动仍较为剧烈.图6和7的结果表明:时均速度较时均压力对湍流模型的改变更敏感,尤其是近壁面流动,如图6(a)、7(a)所示,这是不同湍流模型调用不同壁面解析函数造成的.图6(b)、7(b)中不同湍流模型得到的流场流速和压力的周向分布基本相同,说明在远流场各个湍流模型预报的数值仿真结果基本相同,但图7(b)中 x/L=0.90处的流速和压力分布仍有相对明显的波动,这是由于艇体尾部逆压力梯度造成尾部流速降低,边界层增厚,从而改变流场流速和压力基本对应的变化趋势.

图5

图5 流场提取截面

Fig.5 Slices of flow field

图6

图6 流场时均轴向速度的周向分布

Fig.6 Circumferential distribution of time-averaged axial velocity

图7

图7 C_p的周向分布

Fig.7 Circumferential distribution of C_p

据上推断,10种湍流模型在Sub-off流场时均变量的捕捉上存在差异,且主要应该体现在流场黏性上,进而导致其他流场变量的变化,在本文采用的时间及空间分辨率的情况下,LES-DE、LES-W同比效果更差,原因可能是空间或时间分辨率不足导致的湍流模型对流场黏性特征的预报存在差异.

不同湍流模型对潜艇周围流场参数的平均量预测结果存在差异,同时对流场的二次量如雷诺应力的预测结果亦有不同.图8所示为不同湍流模型仿真得到的潜艇周围流场的雷诺正应力RNS,其余湍流模型的结果由于数值过小,出于文章简洁的考虑未在文中列出.如同文献[17],LES模型对流场二次特征的捕捉较好,DES模型次之,RANS模型最差,故研究流场二次量时模型的选择会引起数值仿真结果的巨大差异.

图8

图8 艇体周围雷诺正应力预报结果

Fig.8 Forecast results of RNS

综上可知,就模型尺度下的潜艇而言,RANS类湍流模型不能有效地捕捉潜艇流场的流动细节,对湍流应力有明显的抑制;DES类湍流模型普遍会抑制潜艇尾流的流动涡结构,同时预报的宏观力亦存在明显的偏差;LES类湍流模型在捕捉潜艇流动涡结构、湍流应力方面表现优异,但是除LES-S模型外,针对宏观力的预测均出现明显的误差.因此在开展潜艇水动力及流场特征分析的综合研究中,LES-S模型是一个较为优异的选择.

4.3 潜艇周围涡场特征

基于上述分析得到的数值仿真策略对Sub-off的流场特征进行计算.由于艇体表面附体的存在,使得其周围流场较为复杂.图9所示为潜艇周围的涡系简图,涡系通过Q标准的重新采样体积进行表达,Q标准是Hunt等^[19]提出的一项涡的判断方法.其中附着在围壳或尾舵附近的稍涡、马蹄涡(项链涡)会改变螺旋桨的来流均匀性而增加螺旋桨的载荷波动、提高空化风险及增加水动力噪声,湍流度巨大的涡结构会提高艇体表面激励力,引发艇身结构振动,削弱潜艇的隐身性能.艇体附着涡是由边界层流动引起的速度梯度变化形成,在边界层较为稳定的平行中体段强度微弱,在几何特征变化明显的首尾段有较强的艇体附着涡.图10所示为围壳附近多截面流场变量,图中:w_i,w_j,w_k为涡量在3个方向的分量,图11对比围壳和尾舵表面的压力变化.从图10可以看出,围壳在轴向首先出现正压力梯度, 在尾段出现负压力梯度,而其附近的涡量表明稍涡的脱落和压力梯度的变化存在密切关系,同时值得注意的是w_k可能是稍涡变化的重要变量,存在明显脱落和分离的w_i、w_j是围壳随边涡变化产生的重要变量.由于尾舵和围壳相似的几何特征,本文认为尾舵稍涡和围壳稍涡的形成及影响因素相似,但尾舵稍涡的强度明显大于围壳稍涡.从图11可以看出,围壳和尾舵顶面和侧面的压力梯度差异可能是导致强度不同的关键因素, 需要注意的是图10和11中均对变量进行无因次化处理以方便对比分析.文献[17]的研究表明马蹄涡的存在与围壳(尾舵)、艇体边界层的三维流动有关,本文不对其进行深入研究.

图9

图9 Sub-off周围涡系

Fig.9 Vortex-system around Sub-off

图10

图10 围壳附近多截面流场变量

Fig.10 Multi-slice flow field variables near the sail

图11

图11 围壳和尾舵表面C_p分布

Fig.11 Distributions of C_p on surface of sail and fin

5 结语

本文以Sub-off为几何模型,针对特定工况,就10种湍流模型计算获得的结果进行分析,得出如下结论:

通过空间及时间收敛性分析,基于理查森外推法判定本文选用的网格方案具有良好的网格收敛性,最终采用了合适的网格开展后续的数值计算.本文选用的湍流模型中,LES-S模型在潜艇的宏观量如阻力、时均流场、伴流方面具有良好的性能或“随群性”;同时,在流场的二次量预报方面,LES类模型获得的结果较其他模型更加丰富.

潜艇涡系因其表面附体的存在较为复杂,但由于围壳和尾舵的几何相似性,其涡系具有内部形成机理的相似性,但因为附体装配位置及参数的差异,尾舵的涡系耦合和演化有较大差异.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

毕毅, 高霄鹏, 王波, 等.

潜艇水动力噪声的自航模试验技术研究

[J]. 海军工程大学学报, 2007, 19(5): 40-43.

, GAO

Xiaopeng

, WANG

, et al.

Acoustic test technique for submarine using remote submersible model

[J]. Journal of Naval University of Engineering, 2007, 19(5): 40-43.

DOI:10.1016/j.ijnaoe.2019.04.002 URL [本文引用: 1]

[2]

PAN Y

, ZHANG H

, ZHOU Q

Numerical simulation of unsteady propeller force for a submarine in straight ahead sailing and steady diving maneuver

[J]. International Journal of Naval Architecture and Ocean Engineering, 2019, 11(2): 899-913.

[3]

POSA

, BALARAS

A numerical investigation of the wake of an axisymmetric body with appendages

[J]. Journal of Fluid Mechanics, 2016, 792: 470-498.

DOI:10.1017/jfm.2016.47 URL [本文引用: 1]

[4]

YAO H

, ZHANG H

, LIU H

, et al.

Numerical study of flow-excited noise of a submarine with full appendages considering fluid structure interaction using the boundary element method

[J]. Engineering Analysis with Boundary Elements, 2017, 77: 1-9.

DOI:10.1016/j.enganabound.2016.12.012 URL [本文引用: 1]

[5]

王超, 郑小龙, 李亮, 等.

Y +值对潜艇流场大涡模拟计算精度的影响

[J]. 华中科技大学学报(自然科学版), 2015, 43(4): 79-83.

WANG

Chao

, ZHENG

Xiaolong

, LI

Liang

, et al.

Influence of Y + on the calculation of submarine flow field characteristics of LES calculation accuracy

[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2015, 43(4): 79-83.

DOI:10.1016/j.ijheatfluidflow.2007.05.007 URL [本文引用: 1]

[6]

, XIAO Z

, CHEN H

, et al.

Simulation of wing-body junction flows with hybrid RANS/LES methods

[J]. International Journal of Heat and Fluid Flow, 2007, 28(6): 1379-1390.

[7]

李士强, 肖昌润, 曹植珺.

基于STAR-CCM+的潜艇尾流场及水动力数值分析

[J]. 中国舰船研究, 2018, 13(Sup.1): 29-35.

Shiqiang

, XIAO

Changrun

, CAO

Zhijun

Numerical analysis of wake flow and hydrodynamics for a submarine based on STAR-CCM+

[J]. Chinese Journal of Ship Research, 2018, 13(Sup.1): 29-35.

[8]

李士强, 肖昌润.

基于STAR-CCM+的潜艇舵翼水动力性能研究

[J]. 舰船科学技术, 2019, 41(11): 37-42.

[本文引用: 2]

Shiqiang

, XIAO

Changrun

Research on hydrodynamic performance of submarine rudder wing based on STAR-CCM+

[J]. Ship Science and Technology, 2019, 41(11): 37-42.

[本文引用: 2]

[9]

GROVES

, HUANG

, CHANG

Geometric characteristics of DARAPA Sub-off models

[R]. Maryland, USA: David Taylor Naval Ship R&D Center, 1989.

[10]

, THOMAS

Summary of DARPA Sub-off experimental program data

[R]. Maryland, USA: David Taylor Naval Ship R&D Center, 1998.

[11]

HUANG

, LIU

, GROVES

, et al.

Measurements of flows over an axisymmetric body with various appendages in a wind tunnel: The DARPA Sub-off experimental program

[J]. [C]//Proceedings of 19th Symposium on Naval Hydrodynamics. Seoul, Korea: National Academic Press, 1992.

DOI:10.1017/jfm.2017.20 URL [本文引用: 1]

[12]

KUMAR

, MAHSH

Large eddy simulation of propeller wake instabilities

[J]. Journal of Fluid Mechanics, 2017, 814: 361-396.

[13]

胡健, 耿冲, 冯峰.

基于大涡模拟的螺旋桨梢涡数值分析

[J]. 华中科技大学学报(自然科学版), 2017, 45(11): 68-73.

Jian

, GENG

Chong

, FENG

Feng

Numerical analysis of propeller tip vortex based on large eddy simulation

[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2017, 45(11): 68-73.