航行器水面滑行流场特性数值仿真

刘富强,孙 元,王广平,王雪峰

(中国船舶集团有限公司 第705 研究所,陕西 西安,710077)

0 引言

鱼雷作为我国海军水下攻防的主战装备,因其具有水下隐蔽性好、打击威力大等特点,已成为海军的重要打击力量。传统鱼雷由于水下航行阻力大、质量限制等原因,导致其航程较小,无法有效突破敌方防御[1-2]。空中导弹具有阻力小、航程长等优点,然而随着世界各国空中拦截水平的提升,导弹极易被敌方雷达锁定,在未命中目标之前即可被摧毁[3-4]。海面低空区域海况较为复杂,当武器在水面滑行或近水面掠水航行时,不易被敌方雷达发现,极大增强其突防能力,低空近水面区域的武器研究已成为各国研究的热点,因此研究一种水面滑行的攻击武器非常必要[5-6]。水面导弹/鱼雷武器运动问题可被简化为航行器水面滑行问题进行研究。

近年来,国内外学者针对舰船、滑行艇等航行器水面滑行问题开展了多项研究。Moctar 等[7]对杜伊斯堡测试船型(Duisburg test case,DTC)在不同速度下的阻力特性进行试验和数值仿真研究,数值计算结果与试验结果吻合性较好;De Marco等[8]通过数值仿真阶梯形船体在静水面的滑水过程,研究了船体的沾水区域特性以及尾部涡流特性。此外,国内外学者对小尺寸圆柱棒水面滑行问题进行了相关研究。Timothy等[9]利用水池进行拖曳试验,对圆柱棒在多种攻角和不同淹没深度下的滑水问题进行了一系列试验研究,获得流体动力的经验解;张珂等[10]对回转体在水面滑行问题进行数值仿真,研究了滑行过程中回转体沾湿区域的流场特性以及空泡特性;张新彬[11]将机器人水面运动简化为细长圆柱棒滑水问题,建立了细长圆柱棒与水面接触相互作用分析模型,研究其运动过程中的力学特性;刘富强等[12]通过对回转体在不同速度、不同淹没深度下水面自由航行过程进行数值仿真,探究其水面滑行流场特性和流体动力特性。

上述研究主要是针对低速航行器,而对于大尺度航行器高速水面滑行问题的研究较少。对大尺度航行器水面滑行问题进行试验研究较为困难,因此提出一种数值仿真方法研究航行器水面滑行尤为必要。文中基于STAR-CCM+数值仿真软件,选用剪切应力传输(shear stress transport,SST)k-ω湍流模型,采用流体体积(volume of fluid,VOF)波和多重参考系(multiple reference frame,MRF)模型运动参考系,构建了航行器水面滑行数值仿真模型。对航行器在不同速度水面滑行工况进行数值仿真,研究其滑行过程中的流场特性和流体动力特性。

1 数值仿真模型

1.1 数值方法

水面滑行问题属于气液两相问题,涉及两相之间的相互作用,在湍流的非直接数值仿真中,应用最广泛的是时均性质的RANS 方程,控制方程包括连续性方程、动量方程、能量方程和相体积方程[13-14]。湍流模型采用SSTk-ω模型[15],其基于Baselinek-ω(BSL)湍流模型,考虑了湍流剪切力的输运问题。采用Schnerr &Sauer[16]空化模型仿真航行器在水面滑行过程中的空化现象。采用MRF运动参考系技术[14,17-18]将航行器运动规律以平移速度的方式赋予参考系,然后把相对速度代入控制方程进行计算。对于滑水问题的求解,为了获得更好的流场特性以及流体动力特性,采用VOF方法仿真气液两相界面问题,可以更好的观察自由液面的变化[19]。

1.2 SST k-ω 湍流模型

SSTk-ω湍流模型基于Baselinek-ω湍流模型,考虑了湍流剪切力的输运问题,通过对涡粘性方程设置限制器,对逆压梯度条件下流动分离的发生和发展具有较高精度的预测。Menter[15]认为BSL 模型的缺陷主要是由对涡粘性的预测不当引起的,可通过给涡粘性设置限制器的方法解决,并且由此提出了SSTk-ω湍流模型。

1.3 Schnerr &Sauer 空化模型

Schnerr 和Sauer[16]在2001 年提出了 Schnerr &Sauer 空化模型,该模型将汽相的体积分数和单位体积液体所含有的空泡数量联系起来,模型对于相间质量传递的表达式为

1.4 MRF 运动参考系

MRF 模型是学者Issa 等[18]于1994 年提出的一种对于多区域计算较为简单的定常计算模型。定义坐标系ox′y′z′固连于绕流物体,并且相对于惯性坐标系oxeyeze(以地面坐标系为例),坐标系ox′y′z′的原点位置矢量为r0,平移速度为vt,角速度为 ω。对于计算域内的任一流体微元,假设其相对于坐标系ox′y′z′原点位置为r,则该点在惯性坐标系oxeyeze中的速度可表示为

式中:v为 流体微元在惯性参考系oxeyeze中的绝对速度;vr为流体微元在运动参考系ox′y′z′的相对速度。

采用运动参考系方法,将绕流物体的运动规律以平移速度或旋转速度的方式赋予参考系,然后把相对速度代入控制方程进行计算。相对速度形式的动量守恒方程需要添加额外的体积力,主要用于描述旋转和加速等非惯性运动,包括科氏加速度、向心加速度等。运动参考系中流体质点受到的体积力为

式中:α和a分别为运动参考系相对于惯性参考系的旋转加速度和平移加速度;等号右边5 个加数项依次为科氏加速度、向心加速度、角加速度、平移加速度和重力加速度。

1.5 VOF 波

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

2 可行性验证

2.1 回转体水面滑行

Timothy 等[9]利用水池进行拖曳试验,对直径为88.9 mm 和165.1 mm 圆柱棒的滑水问题进行了一系列试验研究。基于运动体静水面滑行的数值仿真模型和文献[9]中直径88.9 mm、长度142 2 mm的圆柱棒模型,对圆柱棒以12.19 m/s 水平速度、6°攻角水面滑行工况进行数值仿真。

图1 展示了圆柱棒在沾湿长度为4 倍直径滑水状态下的试验照片和数值仿真计算结果。观察仿真结果发现,在滑水过程中,运动体尾部飞溅产生大量水花,与试验照片吻合性较好。在数值计算过程中,对圆柱棒滑水的升阻力进行监测,该工况下数值计算升力为19.37 N,阻力为10.21 N。

图1 滑水试验照片与数值仿真结果Fig.1 Planing test photo and numerical simulation result

基于式(5)[9]对升阻力特性进行无量纲化表示,升力系数CL为0.032 9,阻力系数CD为0.017 3。对比文献中得到的升力系数0.033 8,阻力系数0.017,升力系数误差为−2.5%,阻力系数误差为2.28%,均小于5%,满足工程误差要求。

式中:FL和FD分别表示运动体航行过程中受到的升力和阻力;ρ表示运动环境介质的密度,文中取水的密度ρ=998 kg/m3;v和d分别为运动体的运动速度和直径。

2.2 航行器水下航行

文中所研究的航行器是以某533 mm 口径鱼雷为原型的去鳍舵简化模型,下文对该简化模型在水下20 m/s 零攻角航行工况进行数值仿真。图2为鱼雷简化模型在水下航行的流场图,观察压力云图可知,鱼雷头部压力达到最高0.3 MPa,全域最低压力为19 182 Pa,明显高于水的饱和蒸气压3 170 Pa,在该工况下不发生空化;由图2(b)速度云图可见,鱼雷头部形成驻点,速度为零,流域最高速度为22.99 m/s,为航行器航行速度的1.15 倍。鱼雷为回转体模型,在不考虑重力情况下流场上下对称分布。

图2 航行器20 m/s 时水下航行流场Fig.2 Underwater flow field around the vehicle navigating at 20 m/s

在航行器水下航行数值仿真计算过程中,对航行器航行阻力进行监测,航行器以20 m/s 在水下稳态航行时的半域阻力为2 372 N,则全域阻力为4 744 N,根据式(6)计算阻力系数为0.106 3。由文献[20]查询得该鱼雷X 型鳍舵布局阻力系数为0.141,其中单一鳍舵阻力系数为0.009,则鱼雷光体阻力系数约为0.105,数值仿真计算阻力系数与文献值误差为1.24%,小于5%,处于工程应用误差范围内,该数值计算结果可信。

其中,A为航行器最大横截面积。

参考经典文献中运动体在静水面滑行工况以及某533 mm 口径鱼雷水下航行工况,对其进行相同工况的数值仿真计算,将数值仿真计算结果与文献中结果相对比,发现数值仿真流场与试验照片吻合性较好,且监测流体动力数值与试验值误差小于5%,在工程误差范围内。因此文中所提出的基于STAR-CCM+对航行器水面滑水航行问题的数值仿真计算方法可行,可用于后续对航行器水面滑水航行问题的研究。

3 数值仿真计算

基于前文提出的航行器水面滑行数值计算模型对航行器在不同速度工况下的滑水特性进行数值仿真,研究航行器在不同速度水面滑行的流场特性和流体动力特性,仿真速度分别为10,20,30,40,50,60,70,80 和90 m/s 共9 组工况。

3.1 计算物理模型

文中对航行器在不同速度水面滑行工况进行数值仿真,其运动体与计算域关于航行器纵平面对称,因此在计算过程中采用半域计算,以提高计算效率。图3 为计算域尺寸以及边界条件示意图,航行器直径D为533 mm,滑水航行攻角取11.5°,航行器在水面滑行初始状态,水面恰好淹没航行器尾端,特征沾湿长度L为2.54 m,v指示航行器水面滑行方向水平向左。

图3 计算域尺寸及边界示意图Fig.3 Diagram of calculation domain size and boundaries

考虑到航行器在水面滑行数值仿真计算过程中流场的充分发展,航行器前端流域长度取2.5L,后端流域长度取5L;高度方向,水平面以上空气部分和水平面以下均取10D;流域宽度方向设置为10D,由于采用对称半域计算,实际流域宽度为20D。

采用STAR-CCM+软件进行数值仿真计算,需要对流域边界进行设定,文中所构建的航行器水面滑行数值计算模型入口采用速度入口,出口采用压力出口,流域上下界以及宽度方向非对称面边界设定为压力出口,对称面边界设定为对称平面,运动体即航行器表面设定为壁面。

在利用计算机软件进行仿真过程中,一般借助网格进行求解。数值计算模型网格的数量和质量会严重影响仿真计算的效率和结果准确性。文中采用切割体网格和棱柱层网格对流域网格进行划分,采用表面重构对网格与结构体衔接处进行优化。同时针对航行器尾端沾水重点区域采用体积控制法进行局部体加密。

在网格模型和计算域保持不变的条件下,采用体形状内不同网格加密尺度对航行器尾部淹没区域进行局部加密,得到网格单元总数目分别为237 万,344 万,496 万,581 万和674 万的5 种计算域网格结果。通过对航行器90 m/s 速度、11.5°攻角、水面恰好淹没航行器尾端特定工况下的数值仿真结果进行对比,验证网格数量无关性。

表1 为不同网格数量航行器特定工况下水面滑行流体动力特性,流体动力系数计算参考式(6),237 万和344 万网格数量模型阻力系数和升力系数计算值明显高于其他网格数量模型计算值;581 万和674 万网格数量模型流体动力系数计算值几乎相同。考虑网格数量较大时,会增加计算机运行负荷,计算耗时增长。因此选用581 万网格数量模型用于后续航行器水面滑行问题数值仿真计算。图4 为581 万网格数量模型航行器尾端网格局部放大图,航行器尾端网格平均尺寸为5 mm,网格质量较好。

表1 不同网格数量下水面滑行流体动力特性Table 1 Hydrodynamic characteristics of the vehicle planing under different mesh numbers

图4 航行器尾部网格局部Fig.4 Local mesh at the tail of the vehicle

3.2 水面滑行流场特性

图5 为航行器不同速度滑行工况水面气液交界面及空化效果图,图中黄色为航行器本体,绿色为航行器在水面滑行时体积分数为0.5 的液面,用于表示气液交界面,红色椭圆框内为尾部空化区域。图6 为航行器在不同速度水面滑行工况下的密度云图,观察发现航行器在水面滑行过程中尾部发生沾湿,兴起波浪,液体飞溅形成水花,并且随着航行速度的提高,兴起水花长度明显增长,波浪浪高明显增高。

图5 不同速度下航行器滑行气液交界面及空化效果图Fig.5 Diagram of gas-liquid interfaces and cavitation effects of the vehicle planning at different speeds

图6 不同速度下航行器水面滑行密度云图Fig.6 Density contours of the vehicle planning at different speeds

观察航行器在滑水速度为10 和20 m/s 时水面效果图(图5)发现,航行器尾部完全被液面包围,几乎不发生空化,无明显空泡;当滑行速度高于30 m/s时航行器尾部出现空化泡,并且随着航行速度的进一步增大,航行器尾部空化区域逐渐增大。对比图6 水面滑行密度云图,在10 m/s 水面滑行速度时,航行器尾部几乎没有空化发生;在50 m/s 水面滑行速度时,航行器尾部有明显的空泡发生,且空泡闭合;在90 m/s 水面滑行速度时,航行器尾部空泡与大气贯通,空泡发生破裂。对比航行器不同速度水面滑行尾部液体飞溅情况发现,随着航行速度的增加,液滴飞溅越高,从最初的小液滴演变为大的尾流波浪。

图7 为航行器在不同速度水面滑行工况下对称面的压力云图分布,观察不同速度下水面滑行最高压力发现:航行器在水面滑水航行过程中,压力最高点位于航行器沾水最前端,此处由于水面滑行冲击作用形成高压区,航行器滑行速度越高,局部最高压力越大。由图可见,航行器水面滑行过程中的低压区主要存在于2 个位置,其一为航行器尾部空化区域,航行器尾部速度略高于航行速度,此处压力较低,易发生空化;其二为在航行器贴近尾部尾流区域形成的锥状低压区。在航行器水面滑行工况数值仿真计算过程中,水的饱和蒸气压为3 170 Pa,航行器10 m/s 工况压力场最低压力为99 186 Pa,压力远高于水的饱和蒸气压,不发生空化;而30 m/s 航行器滑水工况压力场最低压力为3 169.6 Pa,航行器尾部开始空化,与水面滑行密度场以及空化效果图显示结果一致。

图7 不同速度下航行器水面滑行压力云图Fig.7 Pressure contours of the vehicle planing at different speeds

然而观察航行器在以速度为40 和50 m/s 工况下空化效果图发现,航行器贴近尾部尾流区域有锥状封闭区域,且观察50 m/s 水面滑行压力云图发现在该处压力明显小于水的饱和蒸气压。图8 为航行器以50 m/s 水面滑行过程中尾部流线局部放大图,对照该速度下的水面效果、压力场及流线图发现,在航行器尾部出现锥状低压区,低压区内主要包含空化的水蒸气,其密闭于气泡内未与大气相通,泡内出现绕流现象。当航行器水面滑行速度继续提高,在高速滑行工况中,尾流呈波浪状,原低压区与空气相通,空泡溃灭,低压区消失。观察航行器在70 和90 m/s 水面滑行工况的压力场发现最低压力低于水的饱和蒸气压,并且水面滑行速度越高,低压区范围越大,航行器尾部空化越明显。

图8 速度50 m/s 时航行器水面滑行尾部流线Fig.8 Streamline diagram of the vehicle planing at 50 m/s

3.3 水面滑行流体动力特性

在对航行器不同速度下水面滑行过程进行数值仿真的过程中,对其稳定状态下水面滑行流体动力进行监测,监测数据如表2所示。

表2 不同速度下航行器水面滑行流体动力特性Table 2 Hydrodynamic characteristics of the vehicle planing at different speeds

图9 为航行器在不同速度下水面滑行流体动力系数曲线,对比发现航行器在不同速度下的流体动力特性明显不同。当航行器水面滑行速度为20,30,40 和50 m/s 时,阻力系数明显较高,升力系数为负值。其他滑行速度工况下,航行器阻力系数基本相同,升力系数为正值。

图9 航行器水面滑行流体动力特性曲线Fig.9 Curves of hydrodynamic characteristics of vehicle planing

由图5～图7 可见,当滑行速度为20,30,40 和50 m/s 时,航行器尾部由于液滴飞溅,沾湿面积明显大于其他滑行速度工况。图10 为航行器在水面滑行速度为10,50 和90 m/s 工况下的沾湿效果图。观察发现,航行器在水面滑行速度为50 m/s时,表面沾湿区域明显大于10 和90 m/s 工况下的表面沾湿面积。图11 为航行器在水面滑行过程中受到的压力和摩擦力在垂直和水平方向的力分解示意图,水平方向合力为水面滑行阻力,垂直方向合力为水面滑行升力。从图中可以看出,在航行器水面滑行过程中,航行阻力主要由摩擦阻力构成,航行升力主要由压差升力构成。因此航行器在20～50 m/s 速度区间内水面滑行阻力系数明显较高。

图10 航行器水面滑行沾湿效果图Fig.10 Diagram of wetting effect of vehicle planing

图11 航行器水面滑行力分解示意图Fig.11 Diagram of force resolution for vehicle planing

当航行器水面滑行速度为20,30,40 和50 m/s时,升力系数为负值,由图5 可见,此时航行器尾部空泡逐渐生成。图12 为航行器在50 m/s 速度下水面滑行过程中距离尾部1 倍直径处横截面密度云图,可以明显看到在航行器尾部周围有明显的空泡存在,空泡几乎完全包裹着航行器尾部。航行器在水面滑行的过程中,尾端发生空化,形成闭合空泡,空泡内压力接近水的饱和压力。而航行器尾端上部处于沾湿状态或裸露状态,其表面压力接近于大气压。水的饱和压力远小于大气压力,约为大气压力的1/30。因此,在航行器尾端压差作用效果为向下的负升力。

图12 速度50 m/s 时航行器水面滑行截面密度云图Fig.12 Cross-sectional density contour of the vehicle planning at 50 m/s

图13 为航行器水面滑行速度为30,40 和50 m/s 时表面压力云图,航行器表面最低压力点位于尾端下部,约为水的饱和压力,而航行器尾端上方压力为大气压。随着航行器滑行速度的提高,航行器尾端下部低压区面积逐渐增大,航行器尾端负升力显著增强。

图13 不同速度下航行器表面压力云图Fig.13 Surface pressure contours of the vehicle at different speeds

然而随着航行器水面滑行速度的提高,空泡溃灭,与大气相通,泡内压力恢复至大气压,航行器尾部负升力消失。由航行器在70,80 和90 m/s 水面滑行工况下的升力系数可知,随着速度的提高,升力系数略微增大,这与航行器在水面高速航行过程中尾部上端液滴飞溅减少有关。

4 结束语

文中基于STAR-CCM+数值仿真软件,构建了航行器水面滑行的数值仿真模型,对航行器在不同速度下的水面滑行工况进行仿真,通过研究航行器水面滑行过程中的流场特性和流体动力特性,获得一种研究航行器水面滑行流场特性的预报方法,主要得到以下结论。

1) 航行器在水面滑行过程中,低速工况下几乎不发生空化,当速度高于30 m/s 时,在航行器尾端发生空化。此时空泡内压力低于航行器尾端沾湿面压力,空泡发生形变,液面向航行器尾部卷曲,形成飞溅。

2) 航行器在低速水面滑行初始空化的过程中,尾端形成封闭空泡,泡内为低压区,泡内出现绕流,此时航行器的升力为负值,待滑行速度提高,空泡溃灭与大气连通,升力值明显提高。

3) 航行器以不同速度水面滑行时的流场明显不同,导致不同速度下升力系数与阻力系数差异较大,升力甚至出现负值,这主要与不同速度下航行器尾端空化效果不同造成沾湿和表面压力分布存在差异有关。

文中对航行器水面航行流场特性研究采用简化航行器定常研究,在后续研究中将会针对带附体复杂航行器水面动态航行过程进行更深层次研究,包括航行器的跨介质水面航行、多自由度运动特性等,以期服务于航行器水面滑行工程应用。