跨介质是近年来提出的一种新概念,用于航行器、飞行器和潜射类兵器,其主要特点在于可以跨越界面从而实现在不同介质环境中的有效运动,能够隐蔽高效突防、执行多重任务^[1]。此概念可追溯到20世纪中期苏联、美国先后都曾开展过的“飞行潜艇”项目,据报道,2008年,美国国防高级研究计划局(DARPA)重启此项计划,这一新式军事装置将能像潜艇一样在水中行进,又能像飞机一样在空中飞行,前景十分广阔,但技术上面临着诸如外形、结构、动力和材料等诸多难题,可检索到的文献不多,主要以探索性研究、弹道仿真为主^{[2, 3, 4]}。出水阶段和入水阶段是物体跨介质运动的关键阶段,运动过程中,物体受自由面流动、介质变化及耦合作用等因素影响^{[5, 6]},受力情况复杂,附加质量法虽然是工程上解决此类问题的成熟办法,在水下机器人控制、结构体参数辨识等方面应用较多^{[7, 8]},但跨介质过程中附加质量的时变特性极大地增加了研究的难度,成为其模型构建和运动特性分析的关键制约因素。传统的理论推导法只适用于无穷域简单外形物体^[9]的附加质量求解,水池试验方法周期长、花费高,而数值方法无比尺效应、流场信息丰富,国内外对于自由面流动的数值模拟已经做了大量的研究工作^{[10, 11, 12]},这也为时变附加质量的求解奠定了基础。目前对于运动物体的附加质量研究主要针对水下运动体和船舶类物体^{[13, 14, 15, 16, 17]},而对于跨介质运动物体附加质量的研究比较少,主要包括强迫振荡法及其数值化方法,不足之处在于振荡所引入的外界能量会导致自由液面以及流场的变化,与实际运动环境不符,增大了计算误差。本文针对运动物体跨越自由界面过程中的时变附加质量求解问题,采用约束数值试验手段,提出一种快速有效的计算策略,并以球头柱体为对象研究了其跨介质过程中的附加质量变化规律,可以为新概念飞行器/武器的总体设计及优化提供技术支撑。

以物体水空跨越为例,建立如图 1所示坐标系。

图 1 倾斜出水过程的坐标系 Fig. 1 Coordinate system of aslant separation

图选项

参考文献[7],给出物体倾斜水空跨越动力学方程:

令v_y=0,ω_z=0,物体沿x轴方向运动,方程组可简化为

令v_x=0,ω_z=0,物体沿y轴方向运动,方程组可简化为

令v_x=0,v_y=0,物体绕z轴转动,方程组可简化为

式(4)~式(6)右端后4项代表物体所受的水作用力/力矩,改变加速度与角加速度的方向,单独写出,如式(7)～式(12)所示：

于是可得

由流体力学的知识易知,若能获取物体表面的流体信息,那么通过对物体表面的压强进行面积分即可得出F_x+、F_x-、F_y+和F_y-,通过对物体表面的压强同时进行面积分和距离积分即可得出M_y+、M_y-、M_ω+和M_ω-,附加质量计算公式虽然由出水动力学方程引出,但从过程来看,只要能确定不同运动参数下的流场环境,物体是出水状态还是入水状态均不影响结论的成立,故式(13)～式(16)同样适合于入水过程。

综上可得出附加质量的计算策略:建立运动物体的流场数值模型,设定2组速度相同、加速度相反的运动规律(v,v′,ω,ω′)与(v,-v′,ω,-ω′)分别计算,之后对物体表面的压强进行积分,利用式(13)~式(16)即可得出。

本文的数值模拟包含水和空气2种流体,流动速度远远小于声速,所以认为水和空气均为不可压缩流体,把流动当作两相流,自由面为水与空气的交界面,整个过程是一个非定常过程,为了能够精确地跟踪描述水空自由面,采用应用广泛的流体体积(Volume of Fluid,VOF)方法处理。

含相界面的气-液两相流可用如式(17)和式(18)所示的Navier-Stokes方程统一描述:

定义一个流体体积函数F,使F的值等于一个单元内流体体积与该单元体积之比。如果F=1,说明该单元内充满流体;如果F=0,说明该单元内不含流体;如果0<F<1,说明该单元内含有自由界面。F由如下守恒方程决定:

约束方程为

考虑水的黏性,利用稳定性较好的湍流模型“标准k-ε模型”对控制方程进行封闭^{[18, 19]},其模型方程如下。

湍流动能方程k:

扩散方程ε:

数值计算中,整个计算域由网格组成,物体的运动意味着物面边界的移动,即运动条件需要定义在边界条件上,使用自定义函数(User-Defined Function,UDF)对其进行描述,计算过程中被求解器动态地加载,本文中设定的运动规律(v,v′,ω,ω′)即通过UDF实现。

随着边界的运动,部分网格被不断压缩,部分网格被不断拉伸,压缩和拉伸产生的畸变会影响计算的精度,畸变达到一定程度还会产生负体积网格而导致计算失败。物体的跨介质运动过程在数值环境中就是边界的运动,所以计算过程中需要用到动网格方法,包括smooth、layering和remesh等多种网格更新方法^[20],下一时间步的运动情况是由当前时间步的计算结果确定的,各种时间步的体网格的更新是基于边界条件新的位置。

smooth根据边界点上的已知位移来光滑调整流域内节点的位置,不改变网格的拓扑结构,能保证网格质量,适用于小位移问题,对于大位移和强切变问题仅靠节点松弛不能保证网格质量,而且有可能出现网格相交的情况,网格的品质会发生改变,需要用到remesh,通过设置网格控制参数如网格的最小、最大长度,网格单元的偏斜度、尺度函数参数可以控制网格移动后的品质。

考虑到无穷域球体的附加质量存在理论解析解,也有比较可靠的试验资料,故选用球体作为参考外形。如图 2所示,整体计算域为20 m×10 m×10 m的长方体,球体直径D=1 m,非结构网格能较好地适应球体曲面的变化,故模型的面网格采用三角网格,计算域采用四面体、六面体和楔形体;为了兼顾模型曲面网格的精细度与总体计算网格数量,将整个计算域划分为静止区域和变化区域,划分网格时在质量优先的准则下对静止区域使用六面体、楔形体网格,对变化区域使用四面体网格,在计算过程中,静止区域的网格保持不变,变化区域的网格根据边界的运动进行实时更新,模型壁面网格尺寸为0.03 m,外围边界网格尺寸为0.3 m,最终网格数量149万,网格质量绝大部分超过0.7,网格质量及局部网格如图 3、图 4所示。

图 2 计算域示意图 Fig. 2 Schematic diagram of computational region

图选项

图 3 网格质量分布 Fig. 3 Mesh quality distribution

图选项

图 4 局部网格示意图 Fig. 4 Schematic diagram of regional mesh

图选项

选用压力-速度耦合求解器和couple算法,计算域上表面为压力出口,其他表面为壁面,时间步长Δt=0.001 s,计算无穷域球体附加质量系数时浸深位置为自由面y=2 m处。4组运动规律下的求解结果如表 1所示,(v,v′,-v′)代表球体分别以(v,v′)和(v,-v′)进行运动。

从表 1中可以看出,y轴方向与理论值的最大相对误差为0.700 02%,z轴方向最大相对误差为0.543 39%,均与理论值吻合较好,说明选定的速度与加速度不影响计算结果的正确性。

选择运动参数(v,v′,-v′)=(1,10,-10),计算球体垂直出水过程中不同浸深时y轴方向、z轴方向的附加质量系数,结果如图 5所示。从图 5中可以发现,球体顶端在接触自由面时,y、z轴方向附加质量系数已经分别下降了29.6%和16.2%,符合实际的运动过程,说明方法适合于带自由面的附加质量求解。

图 5 球体附加质量系数曲线 Fig. 5 Added mass coefficient curves of sphere

图选项

物体在跨介质运动过程中,由于出水阶段和入水阶段受自由水面影响,所以附加质量是随时间变化的参量,为研究其变化规律,以某简化球头柱体外形为例,头部为半球(r=0.2 m),中间为圆柱(r=0.2 m,h=1 m),尾部为圆台(r₁=0.2 m,r₂=0.1 m,h=0.3 m),计算域与网格划分结果见图 6和图 7,出水阶段和入水阶段使用同一套网格,在数值计算时只需改变重力加速度的方向即可。

图 6 柱体外形计算域 Fig. 6 Computational region of cylinder shape

图选项

图 7 柱体外形及其局部网格 Fig. 7 Cylinder shape and its regional mesh

图选项

为便于比较,将浸深和附加质量无量纲化,即

选用组合(v,v′,ω,ω′)=(1,5,0.1,1)及(v,-v′,ω,-ω′)=(1,-5,0.1,-1),参照球体附加质量的设置条件进行计算,球头柱体出入水的附加质量系数求解结果如表 2、表 3所示,图 8~图 11为结果对比图。

在本文的计算条件下,忽略物体小位移的影响,进行了6步计算,后续处理可得出2步、3步、4步、5步和6步的平均结果,不同步数平均结果的差异不超过3%,以6步平均结果为代表,直接给出其对比结果。从图 8~图 11中容易看出,无论出水还是入水,1步计算结果与6步计算的平均结果吻合较好,从曲线上看比较一致,原因是选择的外形相对简单,网格质量较高,计算结果比较理想。就实际数据而言,1步计算结果与6步计算的平均结果是有细微差别的,这主要源自于数值模型和多步计算。

图 8 随浸深变化的轴向附加质量系数 Fig. 8 Depth-varying coefficients of axial added mass

图选项

图 9 随浸深变化的附加转动惯量系数 Fig. 9 Depth-varying coefficients of added inertia moment

图选项

图 10 随浸深变化的法向附加质量系数 Fig. 10 Depth-varying coefficients of normal added mass

图选项

图 11 随浸深变化的附加静矩系数 Fig. 11 Depth-varying coefficients of added static moment

图选项

1步计算时物体的浸深反映了真实的自由面,所以严格意义上讲,只有1步计算的结果最能符合物理实际,但1步计算时流场发展不够充分,可能会导致结果偏离真实值,就工程应用而言,还需要考虑外形复杂程度、网格质量以及参数选取等因素对结果的影响。多步计算能够平抑结果的波动,但随着计算步数的增多,物体的浸深位置会发生变化,变化量的大小取决于设定的时间步长和运动速度,从而导致计算误差增大。只要选取适当的计算条件,计算结果就能够反映出实际的运动情况,而且能够降低结果对影响因素的敏感度,所以应用时可以根据情况选择使用。

跨介质的出水和入水过程中,轴向附加质量系数、法向附加质量系数、附加转动惯量系数和附加静距系数随浸深的变化曲线相似,除法向附加质量系数曲线外,其他3个均在浸深0.5附近存在趋缓区域,这与球头柱体的细长外形有关,在出水阶段和入水阶段,浸没于自由水面下的部分随球头柱体运动的规律比较接近,附加质量受此影响从而表现出趋势的一致性。从量级来看,轴向附加质量系数、附加静距系数最大均不超过0.1,附加转动惯量系数最大不超过0.06,相对法向附加质量系数而言比较小,这是因为球头柱体在全浸没状态下时法向的迎水面积远大于轴向的迎水面积,也说明其法向速度在水流体的作用下相对更难改变。

对比出水和入水过程,同样的浸深下,入水时柱体的轴向附加质量系数高于出水,最小超出29%以上,入水时的法向附加质量系数略高于出水时的数据;出水时的附加转动惯量系数与附加静距系数均高于入水时的数据。虽然曲线形状相似,随浸深变化的趋势一致,但数据差别较大,这也从附加质量的角度证明了出水不是入水的逆过程。

本文综合运用动网格、两相流等数值方法,针对跨介质运动物体出水过程和入水过程的时变附加质量问题,提出了一种快速有效的计算策略,推导了附加质量的计算公式并进行了验证分析,以球头柱体外形为例,对其垂直跨越自由面过程中的附加质量随浸深变化的规律进行了研究,结果表明:

1) 新的计算策略能够对跨介质运动物体的时变附加质量进行计算,精度符合工程要求。

2) 跨介质运动物体的出水过程和入水过程的附加质量变化趋势相同,规律相近。

研究结论可应用于跨介质运动体的设计、性能分析预报,对其他带自由面问题如水下发射、空投入水兵器的研究也具有一定的参考价值。限于篇幅,本文未考虑波浪因素,而且运动限定在三自由度的纵向平面内,波浪、多自由度耦合对跨介质运动物体的附加质量产生何种影响是后续的重点研究方向。