超空泡航行体借助超空泡减阻技术, 使其航行阻力大幅降低, 航速得到突破性提高[1]。但由于大范围空泡的包裹使得超空泡航行体流体动力布局形式与水下常规航行体流体动力布局形式存在明显差别[2], 尾舵设计为折叠式的全动舵可满足航行体从全沾湿、局部空泡到超空泡不同航行阶段的操纵需求。
超空泡航行体初始航行时尾舵在航行体内逐步展开为“十字”全动舵, 从而实现航行体运动控制, 此动态过程属于复杂瞬态流动问题, 其流体动力载荷特性关系着尾舵能否正常展开, 直接影响航行体的运动稳定性, 对此有必要深入研究水下航行体尾舵展开流动问题, 为航行体运动稳定性分析以及尾舵展开机构设计提供依据。针对舵/翼展开问题, 航空航天领域具有广泛的研究基础[3-4]。对于水下超空泡技术领域, 佟力永[5]采用数值计算方法研究了单独舵片在不同操舵状态下的非定常流体动力特性。黄闯等[6]构建了研究楔形舵片失速特性的数值计算模型和水洞试验系统, 给出了舵片升阻力特性和低压面空化特性随舵角的变化关系。Xu等[7]通过数值研究手段分析了自由液面状态下超空泡楔型翼在不同深度、不同攻角的空泡形态及流体动力变化规律。从国内外公开文献来看, 舵/翼气动问题的研究手段和研究成果较为丰富, 舵片单自由度操纵的水动力学研究基础也较为广泛, 但是关于舵片复杂动态操纵方面的研究成果则相对较少。基于此发展可行的尾舵动态展开数值计算方法, 实现对水下航行体尾舵展开过程流体动力载荷的研究分析, 具有重要的理论与应用价值。
本文采用基于嵌套网格的数值计算方法, 利用网格域间的相对运动实现对尾舵展开过程的流场求解, 研究了尾舵不同展开速率、不同攻角以及尾舵展开前后的水动力载荷变化规律, 并以可视化手段进行了尾舵受力机理性分析。
1 计算模型及方法 1.1 数值计算方法超空泡航行体尾舵展开过程的数值计算基于有限体积方法, 对雷诺时均Navier-Stokes(RANS)方程进行离散处理, 求解过程涉及到的控制方程为连续性方程和RANS时均方程。湍流脉动对时均流动的作用通过湍流模型来闭环求解, 目前基于k-ω的SST湍流模型, 对于水下流场求解具有适应性强、计算精度高的优势[5], 因此本文计算求解采用SST湍流模型。超空泡航行体尾舵展开过程不涉及自然空化/人工通气空化问题, 故未采用相关多相流模型和空化模型。数值计算涉及的数学方程包括
(1) |
(2) |
(3) |
式中, 湍流黏度μt为k和ω的函数, 即
(4) |
各系数的推荐取值为: σk=2.0, σω=2.0, α=5/9, β=0.075, β*=0.09。
方程中对流项和扩散项均采用二阶迎风格式进行离散, 时间采用二阶离散格式, 时间步长为1×10-3。流场求解采用分离流求解器, 根据SIMPLE算法实现分离流模型的求解更新。
1.2 计算流域及网格超空泡航行体尾舵为“十字”全动舵, 采用楔形超空泡翼型。为实现超空泡航行体尾舵展开过程的运动模拟, 采用网格区域相对运动的思想, 结合重叠网格方法, 通过挖洞和插值方式为各网格区域流场计算提供边界信息的传递服务, 从而建立不同网格区域的耦合关系, 实现整个计算区域的数值求解[8-9]。
基于上述网格策略, 首先建立计算流域, 计算流域整体采用圆柱型, 另外建立尾舵运动区域, 各流域尺寸和边界条件设置为: ①进口:进口面距离航行体头部端面为1L(L为航行体长度), 设置为速度进口, 数值为18 m/s; 圆柱体外径为10D(D为航行体直径), 同样为速度进口边界条件。②出口:出口面距离航行体尾部端面为2L, 设置为压力出口, 数值为1.52×105 Pa。③重叠网格:尾舵运动区域外流域表面, 设置为重叠网格边界, 与航行体外流域联合创建生成重叠网格交界面。④壁面:航行体壁面和尾舵壁面为无滑移壁面条件。
采用切割六面体网格完成整个流域网格剖分, 在航行体以及尾舵壁面进行网格加密, 以捕捉较大梯度变化下的流场信息。近壁面采用混合壁面函数处理方式, 以实现壁面附近的黏流梯度有效求解, 近壁面边界层的第一层网格y+值满足壁面函数要求。两部分区域网格通过重叠网格方式实时构建整体求解网格, 图 1为计算区域网格以及尾舵展开前后重叠后的网格示意图。
2 计算方法验证 2.1 网格独立性验证本文采用数值计算方法对超空泡航行体尾舵展开过程进行仿真计算, 故计算网格对于计算结果的准确性相当重要。为验证网格划分形式的有效性, 本文建立了4种不同疏密程度的计算网格, 网格数量分别为0.52×107, 0.96×107, 1.81×107, 2.33×107, 选取尾舵中间截面压力系数作为网格收敛性的衡量指标。图 2给出了4种网格密度的尾舵中间截面压力系数和摩擦力系数分布。对于本文中的物理模型和计算流域尺度, 4种网格密度条件的尾舵压力系数分布类似, 从尾舵压力系数局部放大图和表面摩擦系数分布来看, Grid3和Grid4压力系数计算结果一致性更好, Grid4网格密度更大。兼顾数值计算的准确性和高效性, 本文选择中等网格密度Grid3, 其网格数量为1.81×107。
2.2 计算策略有效性验证本文数值计算涉及定常和非定常2种状态, 以及基于重叠网格的静止和运动形式, 故有必要针对这几种物理过程进行数值计算方法准确性的验证。常规的椭球体和立方体具有阻力系数理论精确解, 选取长短半轴比为4的椭球体为验证对象, 采用上文提到的网格划分策略和数值计算模型, 分别对其定常、非定常的阻力系数进行求解, 其中非定常阻力系数求解采用相对运动思想, 使内流域沿轴向以一定速度平动, 通过重叠网格方式实现内外流场信息传递, 可求解获得轴向运动的阻力系数, 即非定常运动的阻力系数。本文中尾舵展开过程旋转运动求解与这里的平移运动求解类似, 故本文所选取的验证方法有效。数值计算值与理论值对比结果如表 1所示。计算对比结果可知, 本文建立的重叠网格策略和运动求解模型具有较高的计算精度, 可用于模拟本文所研究的尾舵展开运动过程。
3 计算结果与分析 3.1 尾舵不同展开速率状态水动力载荷分析尾舵展开速率不同, 流体动力响应结果也会有一定差异, 对此本节分析了尾舵90°/s, 180°/s, 360°/s 3种不同展开速率状态下尾舵的水动力载荷变化情况, 其结果如图 3所示。可以看出, 在初始展开位置和展开终了位置附近, 尾舵流体动力特性变化比较缓慢, 特别是在展开速率较小时此现象较为显著; 在尾舵中间展开位置附近, 尾舵阻尼力和阻尼力矩随时间基本呈线性变化规律, 且尾舵展开速率较大时, 其阻尼力和阻尼力矩变化率越大。值得关注的是, 尾舵不同展开速率对于尾舵完全展开时尾舵的阻尼力和阻尼力矩影响较大, 其具体结果见表 2, 尾舵展开速率越高, 尾舵所受阻力减小, 尾舵阻尼力矩随之减小, 尾舵法向力略有增加。同时开展了尾舵完全展开后非定常流体动力计算, 其结果与90°/s展开速率完全展开时非定常流体动力结果接近, 两方面结果均说明了尾舵展开过程非定常影响显著, 尾舵展开速率对于尾舵水动力特性影响较大, 其主要原因在于尾舵展开过程周围流体除轴向运动外还存在随尾舵转动产生的周向运动, 故尾舵展开过程附近流场相对比较复杂, 具体的流场细节分析如下文所述。
展开速率/((°)·s-1) | 阻尼力系数/10-2 | 阻尼力矩系数/10-4 | 法向力系数/10-3 |
360 | 1.08×10 | -2.36×10 | 2.82×10 |
180 | 1.28×10 | -2.59×10 | 2.40×10 |
90 | 1.43×10 | -2.81×10 | 2.25×10 |
瞬态计算 | 1.38×10 | -3.02×10 | 2.18×10 |
图 4为不同展开速率尾舵流域中间截面的速度矢量图, 其中舵梢和舵后速度矢量分布差异较大, 尾舵360°/s展开速率时, 来流在舵梢形成绕流, 舵面转动带动舵后流体随体运动, 尾舵展开速率较小时, 如90°/s, 180°/s, 舵后一部分流体随旋转体运动, 另外一部分随主流方向流向尾舵后方, 流动已基本趋于稳定状态。图 5为不同展开速率尾舵表面压力云图, 可以看出, 尾舵展开速率越高, 尾舵前缘高压区域越大, 特别是在舵梢部分。原因在于尾舵不同展开速率带来附近流场改变程度不同, 随着展开速率的增大, 流体对于舵面的相对线速度提升, 且随着转动半径的增大, 舵梢部位相对线速度最高。尾舵以360°/s的速率展开时, 受舵面流动滞止作用, 由舵根至舵梢形成了逐渐扩大的高压区域, 非定常效应较为明显。而以90°/s速率展开时, 尾舵舵根和舵梢前缘压力较为均衡, 可以认为此时尾舵已基本处于非定常稳定流动状态, 其流体动力响应值也较为稳定, 舵根尾缘由于较大的迎流面积形成低压区域, 其反映出来的阻尼力和阻尼力矩较大。故尾舵展开速率越高, 在尾舵展开末程时所需的尾舵驱动力矩越小, 但维持展开状态的力矩应与非定常稳定流动时的流体动力一致。
3.2 尾舵不同攻角状态水动力载荷分析超空泡航行体由于特殊的鸭式流体动力布局, 流体动力作用中心远离航行体质心, 且位于质心之前[2], 在初始阶段航行体会形成一定正攻角, 故尾舵展开过程耦合攻角的变化可能面临更加复杂的水动力载荷环境, 本节以180°/s的展开速率, 计算了0°, 3°, 6°攻角状态下尾舵的水动力载荷变化情况, 其结果如表 3所示。从结果可以看出, 尾舵展开过程舵面升力系数基本呈线性增长方式, 法向力系数变化幅度不大, 来流对尾舵的阻尼作用随攻角增大而减小, 因而可以认为, 初始状态存在的正攻角更有利于尾舵舵面的展开。
攻角/(°) | 阻尼力系数/10-2 | 阻尼力矩系数/10-4 | 法向力系数/10-3 | 升力系数/10-3 |
0 | 1.28 | -2.59 | 2.40 | 0.063 |
3 | 1.23 | -2.52 | 2.12 | 3.52 |
6 | 8.56 | -1.80 | 1.75 | 7.03 |
为进一步分析舵面受力情况, 图 6给出了舵面展开后舵根位置(距舵根0.15L, L为尾舵展向长度) 和舵梢位置(距舵根0.85L)的压力系数分布。对于舵根位置, 随着攻角增大, 压力面压力逐渐增大, 吸力面前缘形成更为显著的负压力峰, 降低了轴向阻尼力, 舵根尾缘位置迎流效果明显, 非定常作用下形成周期性的旋涡流动, 使得尾缘面压力分布稳定波动。对于舵梢位置, 其压力面和吸力面压力分布特征与舵根位置类似, 不同的是, 舵梢因舵宽变小, 加之周向转动作用影响, 舵梢尾缘面的阻滞作用减弱, 尾缘压力分布较为平缓, 吸力面在攻角范围内属于背流区域, 随着攻角的增大, 产生了低于尾缘面的压力, 抵消了尾舵产生的部分阻力。综合舵根和舵梢截面压力系数情况, 随着攻角的增大, 尾舵前后舵面的压力差减小, 故展开后的水动力载荷减小, 与表 3分析结果一致。
3.3 尾舵展开前后水动力载荷分析除尾舵展开过程外, 尾舵张开前后流体动力特征变化显著, 因此有必要分析尾舵张开前后尾舵部分乃至整个航行体的流体动力差异。为减少计算量, 本节以定常计算方式实现尾舵展开前后2个状态的数值模拟计算。表 4量化给出了航行体和尾舵的力及力矩情况, 明显可以看出, 尾舵展开前后阻力相差较大。
尾舵状态 | 航行体阻力系数 | 舵阻尼力系数 | 舵法向力系数 | 舵阻尼力矩系数 |
未展开 | 2.49×10-1 | -2.51×10-3 | -1.07×10-3 | -2.48×10-5 |
舵完全展开 | 2.82×10-1 | 1.38×10-2 | 2.18×10-3 | -3.02×10-4 |
差值比/% | 11.71 | 118.20 | 149.15 | 91.77 |
尾舵未展开时, 完全折放于尾部凹槽内, 来流与凹槽形成滞止区域, 即在尾舵舵梢部分高压区, 如图 7a)所示, 对于单个舵片而言, 形成了正向推力。尾舵展开后外露于航行体, 尾舵后部形成连续的低压区, 尾舵由舵根至舵梢宽度逐渐较小, 相应舵根部分低压区域更为明显, 如图 7b)所示。相比之下, 尾舵展开后迎流条件使其阻力大幅增加, 增幅为118.2%, 其阻尼力的增加使维持舵的阻尼力矩值也相应提升, 展开后航行体的阻力也明显提升。尾舵展开前尾舵前缘形成的高压区域面积大于尾缘高压区域面积, 因而产生向下的法向力, 尾舵展开后在舵根下部形成流动滞止区, 舵片法向方向产生压差, 形成垂直向上的法向力。
4 结论本文采用数值计算方法开展了超空泡航行体尾舵展开过程流体动力载荷计算研究, 获得以下结论:
1) 尾舵展开速率越高, 在尾舵展开末程时所需的尾舵驱动力矩越小, 原因在于展开速率越高非定常效应增强, 带动舵后流体发生随体运动, 而展开速率较小时非定常流动趋于稳定, 舵根尾缘较大的迎流面积形成低压区域, 使其阻尼力和阻尼力矩较大;
2) 尾舵展开时攻角增大, 舵根位置吸力面负压力峰更为明显, 舵梢吸力面产生了低于尾缘面的压力, 二者降低了尾舵展开后的轴向阻尼力, 阻尼力矩也随之减小;
3) 尾舵展开前后, 尾舵附近流动滞止面变化使其阻力、阻尼力矩大幅提升, 压差作用的改变使得垂直于流向的法向力大小和方向都发生变化。
[1] | NOURI N M, MADOLIAT R, YAHANGARDY Y, et al. A study on the effect of fluctuations of the supercavity parameters[J]. Experimengtal Thermal & Fluid Science, 2015, 60: 188-200. |
[2] |
张宇文, 袁绪龙, 邓飞. 超空泡航行体流体动力学[M]. 北京: 国防工业出版社, 2014.
ZHANG Yuwen, YUAN Xulong, DENG Fei. Fluid dynamics of supercavitating underwater vehicles[M]. Beijing: National Defense Industry Press, 2014. (in Chinese) |
[3] |
范国芳, 李强, 操小龙, 等. 折叠翼面展开过程中气动载荷特性数值模拟[J]. 气体物理, 2017, 2(6): 48-53.
FAN Guofang, LI Qiang, CAO Xiaolong, et al. Numerical study on aerodynamic load of a high-speed aircraft folding wing in its deploying procedure[J]. Physics of Gases, 2017, 2(6): 48-53. (in Chinese) |
[4] |
李东, 李斌, 刘仙名. 空空导弹折叠舵展开过程气动载荷分析[J]. 弹箭与制导学报, 2014, 34(3): 154-156.
LI Dong, LI Bin, LIU Xianming. Analysis on aerodynamic loads during deployment of folded tails of an air-to-air missile[J]. Journa1 of Projectiles, Rockets, Missiles and Guidance, 2014, 34(3): 154-156. (in Chinese) |
[5] |
佟力永. 超空泡航行体楔形舵片流体动力学特性数值模拟[J]. 船舶力学, 2015, 19(7): 765-772.
TONG Liyong. Numerical simulation on hydrodynamic characteristics of rudders of supercavitating vehicle on maneuvering condition[J]. Journal of Ship Mechanics, 2015, 19(7): 765-772. (in Chinese) |
[6] |
黄闯, 党建军, 罗凯, 等. 楔形舵片失速特性的数值模拟和水洞试验[J]. 哈尔滨工业大学学报, 2016, 48(10): 125-130.
HUANG Chuang, DANG Jianjun, LUO Kai, et al. Numerical simulation and water tunnel experiment on the stalling characteristics of wedge rudders[J]. Journal of Harbin Institute of Technology, 2016, 48(10): 125-130. (in Chinese) |
[7] | XU Chang, KHOO Boocheong. Numerical investigation on free surface effect on the supercavitating flow over a low aspect ratio wedge-shaped hydrofoil[J]. Journal of Hydrodynamics, 2020, 32(1): 20-30. DOI:10.1007/s42241-020-0003-7 |
[8] | YANG Dongmei, SHAO Fei, LI Chuanglan, et al. Overlapping gird technique for numerical simulation of a fast-cruising catamaran fitted with active T-Foils[J]. Journa1 of Marine Science and Application, 2019, 18(2): 176-184. |
[9] | CHEN Hao, QIAN Ling, MA Zhihua, et al. application of an overset mesh based numerical wave tank for modelling realistic free-surface hydrodynamic problems[J]. Ocean Engineering, 2019, 176(3): 97-117. |
[10] |
张宇文. 鱼雷外形设计[M]. 西安: 西北工业大学出版社, 1998.
ZHANG Yuwen. Configuration design of torpedo[M]. Xi'an: Northwestern Polytechnical University Press, 1998. (in Chinese) |