回转体滑水航行流体动力特性研究

鱼雷具有水下打击威力大、隐蔽性好等特点，是我国海军水下攻防的主战装备，已成为建设我国强大海军的重要力量。然而，传统鱼雷由于水下航行阻力大、平台空间小及质量限制等原因，导致其航程较小，无法有效突破敌方防御圈，发射后我方平台易被敌方锁定，对发射母艇的威胁较大[1-3]。空中导弹具有航程长、阻力小、杀伤力强等优点，然而随着各国对于空中拦截水平的提升，空中导弹易被敌方锁定，在未命中目标之前将被摧毁[4-5]。而海面低空区域海况较为复杂，当导弹在近水面掠水航行或者水面滑行时，不易被敌方雷达发现，极大增强其突防能力，低空近水面区域的研究现阶段已经成为了各国研究的重点，因此研究一种水面滑行的攻击武器非常必要[6-7]。其与水下鱼雷相比阻力明显减小，航程显著提高，相对于空中导弹，其隐蔽性明显增强。水面导弹/鱼雷武器运动问题可简化为回转体在水面的滑水问题进行研究。

近年来，国内外对于运动体在水面的滑水问题主要集中于对于舰船、滑行艇等的研究。Moctar等[8]对DTC船模在不同速度下的阻力特性进行实验和数值仿真研究，数值计算结果与试验结果吻合性较好；张明霞等[9]基于STAR-CCM+对小水线面三体船数值仿真，得到不同的水下潜体形状、侧体位置的阻力变化规律；Marco等[10]通过数值仿真阶梯形船体在静水面滑水过程研究船体的沾湿区域特性以及尾部涡流特性；Amromin[11]通过优化船体底部设计，采用空化的方法为船体减阻提供了解决方案。

国内外对小回转体水面滑行的研究较少。Timothy等[12]利用水池进行拖曳实验，对圆柱棒在多种攻角、不同淹没深度的滑水问题进行了一系列的实验研究，获得流体动力特性的经验解；Waid等[13]对圆柱棒在空化和非空化情况下，不同波形水面滑水工况进行实验研究，得到空化特性和波面对圆柱棒滑水问题的影响规律；Serebryakov等[14]对圆柱体沿一定迎角滑水进行了实验研究，分析了滑水体与型腔表面之间的相互作用过程，研究了深浅度对超空化的影响。

国内对于回转体研究主要集中在对于回转体的入水空泡绕流以及水面冲击特性研究，何春涛等[15]基于VOF(Volume of Fluid)方法求解气、水两相流动的RANS方程，并结合动网格技术，对回转体垂直入水空泡生成过程，空泡壁面运动特性和空泡表面闭合特性进行了研究；孙凯等[16]采用商用CFD软件，仿真了平头回转体模型在亚音速、跨音速以及超音速状态下的入水过程，得到了模型在不同速度下的入水阻力特性和头部直径、液体可压缩性及空气激波对入水过程的影响；张新彬[17]在对仿水黾机器人研究过程中，将机器人水面运动简化为细长圆柱棒滑水问题，建立了细长圆柱棒与水面接触相互作用分析模型，研究其运动过程中的力学特性。

对于大尺度回转体高速滑水问题的实验研究较为困难，因此提出一种数值仿真的方法研究回转体滑水问题尤为必要。本文基于STAR-CCM+数值仿真软件，选用k-w湍流模型，采用VOF波，构建运动体在静水面滑水数值计算模型。对回转体在不同速度、不同淹没深度滑水工况进行数值模拟，研究其滑水过程中的流场特性和流体动力特性。

本文所有计算均采用西北工业大学航海学院高速航行器创新实验室的数值模拟工作站完成，数值模拟计算软件版本为STAR-CCM+ 14.02.012。滑水问题属于气液两相问题涉及两相之间的相互作用, 数值方法控制方程包括连续性方程、动量方程和相体积分数方程[18-19]。在对湍流的非直接数值模拟中，应用最广泛的是RANS雷诺平均方程，湍流模型采用SST k-w模型[20]，其基于BSL k-w湍流模型，考虑了湍流剪切力的输运问题。采用Schnerr & Sauer空化模型[21]模拟回转体水面滑水航行过程中的空化现象。对于滑水问题的求解为了获得更好的流场特性以及流体动力特性，采用VOF方法模拟气液两相界面问题，其可以更好地观察自由液面的变化[22]。

STAR-CCM+软件提供VOF波建模，VOF波模型用于模拟轻流体和重流体之间的交界面上的表面重力波。此模型通常与六自由度运动模型一起用于海洋应用。STAR-CCM+提供以下VOF波模型：平波、一阶波、五阶波、椭圆余弦波、叠加波、不规则波等。VOF波用于多相模拟，在默认状态下，轻流体默认为空气，重流体默认为水。本文所研究的运动体静水面滑水问题，采用平波模型构建气液交界面，通过设置平波水面上的点、水面方向模拟回转体在静水面滑水航行工况。

本文采用学者Schnerr和Sauer[21]在2001年提出Schnerr & Sauer空化模型，该模型将汽相的体积分数和单位体积液体所含有的空泡数量联系起来，模型对于相间质量传递的表达式为：

式中：ṁ+为蒸发率；ṁ-为凝结率；α为气相体积分数；RB为空泡半径；n为单位体积内空泡数量；pv为液体饱和蒸气压。

本文采用隐式非定常瞬态求解模型仿真运动体滑水过程[23-24], 其计算时间步长的设置非常关键, 设置时间步长过大, 影响仿真结果的准确性, 设置时间步长过小, 计算较为耗时。因此, 一个合适的计算时间步长对于数值仿真计算同样重要。

本文研究的运动体滑水问题, 计算时间步长选用关系式[10]如公式(3)所示, 最大时间步长为1×10-4s。在本文数值仿真工况中, 运动体速度最大值为21.69 m/s, 近壁面网格尺度最小值为3×10-3 m, 根据库朗数计算公式(4)得最大库朗数为0.723, 不大于一般系统设定默认值1.0, 选取时间步长合适。

式中：V表示运动体的速度；l表示在速度方向上运动体的特征长度；Δt为时间步长；Δx为网格尺度。

Moctar等[8]对DTC船模在6种不同速度下的阻力特性进行实验, 基于本文提出的数值仿真模型和DTC船模, 对6种不同速度下的船模水面航行问题进行数值模拟。

DTC船模具有严格的对称性, 在滑水计算过程中采用半域计算, 减小计算工作量, 对称面附近网格如图 1所示, 对其船模附近流域进行多层加密, 从而确保仿真计算结果的准确性。数值计算采用k-w湍流模型结合速度入口实现船模运动, 外边界设置为速度入口, 出口设置为压力出口。

本节对文献[8]中DTC船模不同速度实验工况进行数值仿真, 其中DTC船模垂线间长5.976 m, 水线宽度0.859 m, 静水面位于龙骨上方0.244 m处, 排水量为0.827 m3, 速度分别为1.335, 1.401, 1.469, 1.535, 1.602, 1.668 m/s。

图 2展示了不同速度下船模稳态航行下的水面示意图。在DTC船模以不同的速度在水面运动数值仿真过程中, 对其船体受到的阻力进行监测, 监测结果如表 1所示, 表中RT表示文献中试验得到的总阻力, RF表示文献中摩擦阻力计算值, RT-S表示数值仿真过程中船体受到的总阻力, RF-S表示数值仿真过程中船体受到的摩擦阻力。将不同速度船模数值仿真模拟结果的计算值与实验值对比, 总阻力误差不超过4%, 摩擦阻力误差不超过1.5%, 符合工程误差范围内。

Timothy等[12]利用水池进行拖曳实验, 对不同直径回转体的滑水问题进行了一系列的实验研究。基于运动体静水面滑水的数值仿真模型和文献[12]中直径88.9 mm、长度1 422 mm回转体, 对回转体以12.19 m/s速度、6°攻角进行一定淹没深度滑水工况进行数值仿真。

图 3展示了回转体在沾湿长度为4倍直径滑水状态下的实验照片和数值仿真计算结果。观察仿真结果发现在滑水过程中, 运动体尾部飞溅形成大量的水花, 与实验照片吻合性较好。在数值计算过程中, 对回转体滑水的升阻力进行监测, 该工况下水面滑行实验和仿真结果对比如表 2所示, 数值模拟升力为19.37 N, 阻力为10.21 N。基于公式(5)对升阻力特性进行无量纲化表示, 升力系数为0.032 9, 阻力系数为0.017 3。对比文献[12]中得到的半域升力系数0.033 8, 阻力系数0.017 0, 升力系数误差为-2.66%, 阻力系数误差为1.76%, 升阻力系数误差均小于5%, 满足工程误差要求。

式中：Fl和Fd分别表示回转体滑水过程中受到的升力和阻力；ρ表示运动环境介质的密度, 本文中取水的密度998 kg/m3; v和d分别为回转体的运动速度与直径。

本章参考国外运动体在静水面滑水的文献, 对其实验工况进行相同工况的数值仿真计算, 将数值仿真计算结果与文献中结果相对比, 发现数值模拟流场与实验照片吻合性较好, 且监测流体动力与实验值误差小于5%, 在工程误差范围内。因此本文所提出的基于STAR-CCM+软件对运动体水面滑水航行问题数值模拟计算方法可行, 可用于后续对于回转体在水面滑水航行问题的研究。

小运动体低速航行可以依靠实验进行研究, 而对于大尺度、高航速的运动体水面滑行问题实验较为困难。本章基于前文提出的运动体滑水数值计算模型, 对回转体在不同工况下的滑水问题进行深入的数值模拟, 研究回转体在不同速度、不同淹没深度滑水工况下的流场特性和流体动力特性。

本章主要对回转体在水面滑水航行不同工况进行数值仿真计算, 其运动体和计算域关于对称面对称, 因此在计算过程中采用半域计算, 提高计算效率。图 4为计算域尺寸以及边界条件示意图, 其中回转体直径为D=88.9 mm, 回转体滑水攻角τ为6°, 回转体长度为L=1 422 mm, v指示回转体的滑水方向水平向右。考虑滑水运动回转体周围流场的充分发展, 设回转体水面沾湿长度Lw为特征长度, 回转体迎水面沾湿点距离速度入口流域长度为6Lw, 回转体尾部距离压力出口流域长度为10Lw; 高度方向, 水面以上空气流域高度设定为8D, 水面以下水深设定为12D; 流域宽度设定为10D, 由于采用半域数值仿真计算, 其实际流域宽度为20D。

采用STAR-CCM+软件进行数值仿真计算, 需要对流域边界设定边界条件, 本文研究的回转体滑水航行模型入口采用速度入口, 出口采用压力出口, 流域上下界以及宽度方向非对称面边界设定为速度入口, 对称面边界设定为对称平面, 运动体即回转体表面设定为壁面。数值计算过程中, 回转体表面速度设定为0, 通过速度入口给定流域速度实现回转体的相对运动。

在利用计算机软件进行仿真过程中, 一般借助网格进行求解。数值计算模型网格数量和质量严重影响仿真计算的效率和结果准确性。在网格边界层底层网格高度、网格模型、计算域保持不变条件下, 对回转体尾部淹没区域局部加密得到网格单元总数目分别为109万, 172万, 198万, 239万和280万的5种计算域网格结果。通过对比回转体12.19 m/s速度、6°攻角特定工况数值模拟结果, 验证网格无关性。

表 3为不同网格数量回转体特定工况水面滑行流体动力特性, 与文献[12]阻力系数0.017 0、升力系数0.033 8相比, 109万, 172万和198万网格数量模型阻力系数计算误差较大, 超过20%;而239万和280万网格数量模型流体动力系数计算误差明显较小, 均小于3%。考虑网格数量较大时, 增加计算机运行负荷, 计算耗时增长。因此选用239万网格数量模型用于后续回转体滑水问题数值模拟计算。

本节对回转体在不同速度下的滑水工况进行数值仿真模拟, 采用无量纲系数Cv表征速度见公式(6), 采用无量纲系数λ表征回转体的淹没深度见公式(7)。回转体滑水特征速度Cv分别取3, 8, 13, 18和23, 沾湿特征长度λ取4。

式中: v为回转体运动速度; d为回转体直径; g取当地重力加速度。

式中: Lk为回转体在轴向的沾湿长度; d为回转体的直径; h为回转体的淹没深度; τ为回转体的滑水攻角。

本节通过对不同速度下回转体的滑水工况进行数值模拟, 研究其流场特性和流体动力特性。图 5表征回转体在不同速度下滑水水面效果图, 低速滑水时回转体尾部沾水面积较少, 随着滑水速度的提升, 回转体尾部溅起大量水花, 回转体尾部发生大面积沾水, 当Cv大于13时, 明显观察到尾部水花飞溅。图 6为回转体在Cv为8, 13, 23时对称面速度云图, 回转体在滑水过程中, 回转体头部出现驻点; 随即在周围展开出现绕流, 流域速度最大值为回转体滑水速度的1.2倍左右; 回转体尾部在滑水初期局部压力低于水的饱和蒸气压时, 发生空化, 同时当尾部附近液面压力高于尾部局部压力时, 液面发生卷曲, 有波浪兴起, 发生飞溅。

回转体在不同速度下滑水时的流体动力特性如表 4所示, 阻力和升力采用系数化表示, 其中下标f表示摩擦作用下流体动力, p表示压力作用下的流体动力。图 7为回转体水面滑水过程中受到的压力和摩擦力在垂直和水平方向的力分解示意图, 水平方向合力为滑水阻力, 垂直方向合力为滑水升力。对照图 7分析表 4发现, 当Cv≥8时, 滑水阻力中摩擦阻力占总阻力的2/3;滑水升力中压差升力占比超过97%;而Cv为3时, 滑水阻力中压差阻力明显高于摩擦阻力, 升阻力系数明显高于其他速度工况。

观察图 5a), Cv=3回转体滑水水面效果图发现, 在回转体以较小速度滑水时, 滑水面周围液滴飞溅较少, 高于静水面区域回转体几乎不发生沾水, 相比高速滑水, 回转体受到的摩擦力较小; 回转体在低速滑水过程中, 水下完全沾湿部分提供向上的升力, 而在高速滑水过程中由于液滴飞溅在回转体高于静水面的表面受到向下的压力产生负升力, 因此在回转体低速滑水时升力系数明显高于高速滑水时升力系数, 尤其压差作用形成的升力部分。

图 8为回转体在不同速度滑水时的流体动力系数曲线。回转体在Cv≥8时升阻力系数几乎不变, 低速滑水时升力系数和阻力系数明显较高; Cv=3时, 阻力系数高于高速滑水阻力系数20%, 升力系数是高速滑水升力系数的3倍, 其与回转体低速滑水时流场无明显液滴飞溅有关; 回转体高速滑水时, 升阻力特性属于回转体的固有特性, 与滑水速度无关, 一般与回转体滑水航行沾湿特性有关。

本节对不同淹没深度/沾湿长度的回转体滑水工况进行数值仿真, 沾湿特征长度λ取1, 2, 3和4, 回转体滑水速度Cv取13。研究不同淹没深度回转体滑水过程中的流场特性和流体动力特性。

图 9表征了不同沾湿长度工况下回转体滑水过程水面效果图, 发现随着滑水淹没深度的提高, 回转体沾水长度明显增长, 并且在回转体尾部附近水面产生的液滴飞溅效果更加明显。图 10为λ=4工况回转体滑水液滴飞溅图, 其飞溅长度超过1倍的回转体长度。

观察图 11不同沾湿长度回转体滑水流域速度云图发现, 在回转体滑水过程中头部出现驻点, 流场速度最大值位于头部绕流区域。在4种不同沾湿长度回转体滑水工况中, 流场速度最大值几乎相同, 最大值为14.3~14.5 m/s, 为流场速度12.19 m/s的1.17~1.20倍, 其速度场几乎不受沾湿长度的影响。

同时在数值计算仿真的过程中, 对回转体滑水的流体动力进行了监测, 在不同沾湿工况下回转体的升阻力特性如表 5所示。回转体在不同沾湿工况滑水流体动力系数曲线如图 12所示, 其升力系数和阻力系数均随着沾湿长度的增加而增大。观察升阻力系数曲线的线性特性发现, 阻力系数曲线一次线性特性较强, 而升力系数相比较弱。其主要由于回转体在不同沾湿工况滑水过程中, 随着淹没深度和沾湿长度的增加, 回转体周围液滴飞溅增强, 静水面以上部分沾水面积增大, 随即受到的压力作用分布不均匀。而压差作用力分解主要体现为压差升力, 对于阻力贡献较小。因此在不同沾湿工况滑水过程中, 回转体阻力系数随沾湿长度变化线性特性较强, 升力系数随沾湿长度变化线性特性较弱。

本文基于STAR-CCM+数值仿真软件，构建运动体在静水面滑水数值计算模型。用数值仿真的方法对回转体在不同速度和不同淹没深度工况下的滑水过程进行数值仿真，得到研究运动体滑水航行流体动力特性的预报方法，为近水面区域武器的研究提供理论参考，主要得到以下结论：

1) 回转体滑水过程中，速度场分布几乎不受沾湿长度和滑水速度的影响，速度场最大值为滑水速度的1.2倍左右；

2) 回转体在Cv≥8滑水过程中，流体动力系数几乎不变，其不受速度大小的影响；而在Cv=3时，流体动力系数较高，阻力系数高于高速滑水阻力系数20%，升力系数是高速滑水升力系数的3倍，其与回转体低速滑水时流场无明显液滴飞溅有关，造成该现象的根本原因是回转体低速滑水时表面压力和摩擦力作用分布不同；

3) 回转体不同淹没深度滑水过程中，其阻力系数随淹没深度增长线性特性较强，升力系数线性特性较弱，其主要与不同淹没深度下回转体表面沾湿不均从而导致压力分布有关。