超空泡航行体流体动力延迟特性数值模拟研究

刘喜燕, 袁绪龙, 罗凯, 杜兆星

(西北工业大学 航海学院, 陕西 西安 710072)

由于空泡延迟效应[1]的存在,当超空泡航行体改变攻角时,新生成的空泡需经过一段时间才能到达航行体尾部。空泡状态的更新相对于航行体姿态的变化存在一定的时间滞后,沾湿面积改变,从而引起流体动力的变化。空泡与航行体之间相对位置关系的实时变化,导致超空泡航行体机动航行过程中的流体动力展现出极强的复杂性和非定常特性,从而使超空泡航行体弹道预报和稳定控制充满挑战。因此,研究超空泡航行体流体动力的延迟特性,具有十分重要的理论与工程应用价值。

目前国内外学者围绕超空泡航行体机动航行过程中的空泡形态和非定常流体动力特性进行了大量研究。李魁彬等[2]基于空泡截面独立扩张原理、空泡镜像模型和航行体动力学方程,推导了恒定推力作用和匀变速情况下超空泡形态的解析解,探索了记忆效应对空泡形态的影响,获得了空化数、空化器尺寸和加速度对记忆效应的影响规律。周景军等[3-4]基于CFD方法研究了超空泡航行体加速和操舵过程中的流体动力特性,获得了航行体空化器、尾舵操舵、通气量、重力效应以及航行体攻角等参数变化过程中的空泡形态及航行体瞬态流体动力特性变化规律。李杰等[5]通过建立超空泡航行体自由运动的非定常流场求解数值模拟方法,对尾拍流体动力的形成和演化规律进行了研究,发现尾拍力与空泡壁面形状、航行速度及尾部沾湿面积等密切相关。陈诚[6]通过建立考虑空泡延迟效应的超空泡航行器动力学模型,对航行器的尾部位置力特性进行了仿真研究,得到了空化器舵角导致流体动力出现附加攻角的规律。王威、何广华等[7-8]采用数值方法分别研究了垂直方向和水平方向的周期性来流扰动对通气超空泡航行体流体动力特性的影响,分析了周期性来流的波长和波幅对空泡形态和流体动力的影响规律。罗凯等[9]通过建立定量求解通气超空化的三维数值模型,研究了通气量变化对空泡形态的影响,获得了通气量突变会使得空泡特征呈现不同程度时间迟滞特性的规律。

综上所述,公开的文献主要集中在考虑空泡延迟效应对动力学建模和流体动力特性的研究,缺乏对非定常流体动力延迟特性形成机理的研究。目前,课题组[10]已通过搭建连续变攻角测力试验平台,定量分析了航行速度、摆动频率及空化器预置舵角等因素对超空泡航行体连续变攻角过程中延迟时间特性的影响。在此研究基础上,本文通过建立连续变攻角运动数值模型,结合对空泡与流体动力变化的分析,进一步揭示不同影响因素下非定常流体动力延迟特性的形成机理,为超空泡航行体动力学建模和机动航行控制技术提供参考。

1 计算模型与方法

1.1 连续变攻角运动模型建立

图1给出了连续变攻角运动模型坐标系示意图。设置地面坐标系为Oxyz,初始时刻航行体质心位置为地面坐标系原点,且xOy平面与航行体的对称面重合,航行体中轴线与x轴重合,空化器指向x轴负方向。设置航行体的转轴中心Or位于航行体质心O后460 mm处(水洞试验中尾支撑转动位置),旋转轴正方向垂直于对称面向外,建模过程中不考虑重力的影响。

图1 坐标系定义

航行体的摆动规律为

式中:θ为摆动角度;θmax为摆动的最大角度,本文中θmax=2°;f为摆动频率;ωz为角速度。

1.2 控制方程

本文通过数值模拟和水洞试验验证相结合的方法开展研究。受制于试验条件,水洞中的水速最多只能达到12 m/s,形成自然空化能力较弱,因此采用主动通气实现超空泡流动。借鉴目前通气超空泡数值模拟研究经验[11],采用流体体积函数(VOF)多相流模型进行水、汽两相流动问题模拟研究。

1) 连续性方程

(3)

2) 动量方程

(4)

3) 湍流方程

本文选取SSTk-ω湍流模型,它在更广泛的流动领域中具有更好的鲁棒性和求解精度。

式中:k为湍流动能;ω为湍流耗散率;Γk和Γω分别为k和ω的有效扩散系数;Yk和Yω分别为湍流导致的k和ω的耗散;Dω为交叉扩散项;Gk,Gω为湍流产生项。

1.3 计算模型与边界条件设置

航行体模型如图2所示,主要由圆盘空化器、圆锥段、圆柱段以及尾喷管组成。

图2 航行体模型示意图

航行体空化器直径D1=19 mm,航行体圆柱段直径D2=54 mm,尾喷管直径D3=32 mm,全长L=432 mm,锥段长L1=181 mm,航行体质心位置距离空化器前端L2=200 mm。

计算域及边界条件设置如图3所示,计算域为长方体,考虑流域边界效应对空泡发展的影响,结合作者对壁面效应的研究经验[12],计算域尺寸长、宽、高分别为34D2,10D2,20D2,航行体放置于对称面上,空化器头部距离计算域左侧为1Lc(Lc为空泡最大长度),尾喷管端面距离计算域右侧为2Lc。计算域左、上、下面设置为速度入口,前、后面以及航行体设置为壁面,右侧面设置为压力出口。

图3 计算域及边界条件

本文采用HEXPRESS软件进行结构化网格划分,重点对空化器、航行体肩部及尾流空泡区的网格进行加密处理。通过研究3套不同单元数量的网格(252万,110万,42万)对计算精度和效率的影响,发现网格数量由42万增加到110万,最大空泡长度增加5.25%,最大空泡直径增大7.67%,而当网格数量由110万增加到252万,最大空泡长度增加0.53%,最大空泡直径增加0.62%,综合分析,最终确定网格数量为110万,网格划分细节如图4所示。

图4 网格划分细节

1.4 计算方法

本文中超空泡航行体做受迫连续变攻角运动,模型不仅具有相对于水速的向前速度,还同时有转动速度。采用动计算域法[13]实现航行体的位置和姿态更新,用UDF控制超空泡航行体按照给定的正弦规律进行摆动,同时将计算获得的航行体受力结果输出至文件。

超空泡航行体流场与运动耦合仿真基于Fluent软件实现,压力-速度耦合算法选用PISO算法,压力项离散格式选用PRESTO格式,动量、体积分数、湍流动能和湍流耗散率的离散选用一阶迎风格式,计算时间步长保持在2.0×10-4s。

计算流程:0～1 s内,航行体保持静止,使空泡充分发展;1～3 s内航行体连续往复摆动若干周期,然后航行体停止运动;最后0.5 s空泡恢复原状态。

2 计算结果与分析

2.1 连续变攻角测力水洞试验系统

依托西北工业大学高速水洞完成连续变攻角测力试验平台的搭建。试验系统如图5所示,由攻角调节与测量系统、通气流量控制与测量系统、流体动力测量系统、压力测量系统、图像采集系统及模型主体组成。模型以尾支撑方式安装于水洞工作段内的电动转台上,试验开始前,使模型保持在0°攻角位置上。通过LABVIEW编写综合控制软件,由单片机控制步进电机,实现模型按设定规律连续转动。通气管路、测压管路和天平导线从导流罩内孔引出至洞体外,压缩空气经流量控制器接至模型通气管。

图5 试验系统构成

2.2 数值模型验证

1) 试验结果重复性验证

为减小试验中通气流量和水洞水速脉动变化带来的影响,确保结果准确性,试验工况均做2～3组重复试验。图6给出水速8.5 m/s,频率1.5 Hz,预置舵角0°工况时,模型摆动方向受力的数据处理前后对比情况。由图6可以看出,2组数据的重合度较好,试验的重复性良好,可以说明试验结果具有较高的可靠性。

图6 试验重复性结果

2) 数值结果与试验结果对比验证

图7给出了通气流量Qg为80 L/min,空化数为0.05时的试验与仿真计算结果对比。由图7可以看出,参照试验空泡结果,仿真空泡最大直径相对误差为5.3%,空泡最大长度相对误差为7.3%,说明建立的数值模型基本合理。

图7 水洞试验与数值仿真结果对比

图8给出了连续变攻角运动过程中摆动方向所受流体动力相对于运动相位之间的延迟效应示意图。计算各个运动周期T内攻角曲线峰值与流体动力峰值之间的时间差Δt,取多个周期内时间差值Δti求均值作为该工况下的流体动力延迟时间t。为了方便对比,对延迟时间t进行无量纲化处理,定义为延迟系数Cyc,其表达式为

图8 延迟时间示意图

(7)

式中:Tck为无量纲参数;Tck=L/v,其中,L为航行体总长;v为参考速度。

图9给出了来流速度v分别为7.0,7.5,8.0,8.5 m/s、摆动频率f均为1.5 Hz、保持相同通气量80 L/min、预置舵角0°等不同工况水洞试验获得的延迟系数Cyc与数值仿真的对比结果。由图9可以看出,仿真结果与试验结果变化规律一致,延迟系数随着速度增加而增大。数值计算得到的延迟系数与试验结果的最大相对误差为10.4%,进一步证明建立的数值模型具有较好的准确性。

图9 延迟系数Cyc随速度v的变化规律

2.3 速度影响分析

图10给出了水速8.5 m/s,空化数0.05,改变攻角方向,获得攻角范围为-2°～2°,每0.2°一个步长攻角下的航行体升力系数变化曲线。由图10可以看出,在自由攻角以内,航行体所受升力系数始终为零,在攻角αfr=±0.8°,航行体圆柱段处于临界沾湿状态;随着圆柱段沾湿,航行体所受升力系数随攻角增大而增大,当到达临界攻角αcr=±1.7°时,圆锥段开始沾湿并在肩部发展形成附体空泡,使得圆柱段沾湿面积减小,圆锥段沾湿面积开始增大,总沾湿面积继续增大,但升力系数增大的斜率减小。总体而言,升力系数随攻角增加而增大,没有呈现迟滞效应。因此,定常流体动力特性不能用来预报航行体在非定常运动时受到的流体动力。

图10 升力系数随攻角变化曲线

图11给出了水速为8.5 m/s,空化数0.05,摆动频率1.5 Hz,摆动攻角范围-2°～2°,3个摆动周期内航行体升力系数变化曲线。由图11可以看出,随着摆动攻角的周期性变化,航行体的升力系数呈现周期性变化规律。在一个摆动周期内,摆动攻角达到波峰或波谷(如图中b时刻)时,对应此攻角下的升力系数未达到最大值,而在滞后于该时刻出现最大峰值(如图中c时刻),说明航行体在连续变攻角运动过程中,流体动力具有延迟特性。

图11 升力系数随摆动攻角的变化规律

为了分析升力系数曲线的变化细节,提取了航行体连续变攻角运动下第一个周期内典型时刻的空泡特征。图12给出了图11中7个特征时刻下航行体对应的空泡形态。

图12 一个周期内典型时刻空泡形态

由图12可以看出,航行体在a时刻之前上表面处于临界沾湿状态,升力系数呈现出随负向攻角增加而小幅度增大的趋势。a时刻为升力系数转折点,此时对应的航行攻角为-1.88°,结合图10可知,在定攻角航行状态下,航行体自由攻角αfr为±0.8°,说明空泡形态发展相对于攻角改变存在滞后性,这使得流体动力特性呈现出差异性;在b时刻,航行攻角达到最大幅值时(α=-2°),航行体肩部开始出现沾湿,此时升力系数未达到最大值,随后航行体开始回摆,此时迎流面肩部的附体空泡开始形成,当法向力达到极值时,航行体回摆至攻角为α=-1.93°处,见c时刻,升力系数的峰值滞后于攻角的峰值,升力系数呈现迟滞特性,迟滞的相位差为0.043°。在d时刻,航行体升力系数出现拐点,分析认为,这是由于航行体圆柱段大部分被附体空泡包裹,整个航行体尾部沾湿面积的变化率降低;在e时刻,航行体沾湿面积达到最小,见该时刻空泡形态;在t=T/2时刻后,升力系数迅速增加,可知沾湿面积发生了快速的变化,这是因为形成附体空泡向锥段方向发展并与之融合,圆柱段尾端的空泡未能及时发展,空泡发生溃灭,从而使得沾湿面积增大,空泡融合发展如图13所示;在f时刻,航行体圆柱段附体空泡与主空泡发展融合,在闭合过程中产生压力峰,上下表面压力的不对称分布形成了作用于航行体的升力,随后圆柱段上的空泡继续发展,沾湿面逐渐消失,并在g时刻重新包裹整个航行体。

图13 附体空泡融合、发展示意图

表1统计了不同速度下航行体升力系数出现峰/谷值(主峰)和脉冲峰值(副峰)时与攻角的对应关系,其中,Δα为相位角差值,主峰和副峰如图11中标志所示。由表1可以看出,随着速度增大,升力系数主峰值相位均滞后于摆动攻角,且相位角差值增加,升力系数主峰值呈减小趋势,延迟效应加剧;随着速度增加,升力系数副峰值对应的相位角提前,升力系数主峰值与副峰值比值增大,在v=8.5 m/s时,高达2倍。

表1 升力系数与攻角对应关系

2.4 摆动频率影响分析

图14给出了摆动频率分别为f1=1.0 Hz,f2=1.5 Hz,f3=2.0 Hz,速度v=8.5 m/s,空化数σ=0.05条件下,第一个T/2周期内航行体从α=-1°(↓)到向上α=-1°(↑)的运动过程中,相同攻角状态时的空泡形态对比。

由图14可以看出,从α=-1°(↓)摆动到α=-2°(→)过程中,3种摆动频率下的空泡形态发展趋势基本一致,摆动频率的增加使得航行体自由攻角范围增大,如图14所示。α=-1.8°(↓)时的空泡形态,此时仅f1工况的航行体圆柱段出现沾湿;当攻角达到最大幅值攻角(α=-2°)时,沾湿面积有S(f1)>S(f2)>S(f3),且航行体锥柱结合位置开始出现沾湿。随后回摆运动中,航行体圆柱段形成附体空泡,摆动频率的增加加剧了附体空泡的发展,使得航行体沾湿面积减小。攻角由-1.8°(↑)向-1.0°(↑)连续变化过程中,空泡变化呈现明显的迟滞发展特性,并随着频率的增大,迟滞作用减弱。具体体现为:在f1工况下,随着攻角幅值减小,附体空泡在轴向扩张变化率较小,在径向基本未发展变化,使得沾湿面积变化率较小,而其他2种工况下,随着攻角幅值减小,沾湿面积呈现较明显减小。

表2统计了第一个摆动周期内不同摆动频率下升力系数与攻角的对应关系。图15给出了3种摆动频率下升力特性随摆动攻角的变化曲线。结合表2和图15可以看出,在同一摆动频率周期内,航行体升力系数呈周期性变化,随着摆动频率增大,与摆动攻角的相位差Δα减小,升力系数的迟滞特性减小,升力系数的主峰峰值增大。随着摆动频率增加,升力系数副峰在回摆过程中逐步向后移动,其副峰峰值减小。从第一个周期来看,在f1工况下,升力系数副峰峰值出现在航行体负向恢复至零攻角之前,所处攻角为α=-0.53°;在f2工况下,升力系数副峰峰值出现在正向摆动至最大攻角过程中,所处攻角为α=0.19°;在f3工况下,未出现明显的升力系数脉动,但在航行体在正向摆动至最大攻角过程中仍存在一个较小的脉动峰值,所处攻角为α=0.90°。分析认为,摆动频率增加加剧了对空泡流场的扰动[11],使得空泡内外压差发生变化,圆柱段附体空泡发展过程不同,f1和f2工况下圆柱段附体空泡形成了孤立气泡,最终发生溃灭而产生较大的冲击载荷[14],而f3工况的附体空泡发展融合于主空泡,摆动频率的增大促进了空泡的融合。

表2 第一个周期内升力系数与攻角对应关系

图15 升力系数随摆动攻角变化规律

2.5 空化器预置舵角影响分析

图16给出了空化器预置舵角β分别为10°和20°时,航行体0°攻角、空化数为0.05工况下的空泡形态。

图16 航行体0°攻角状态下的空泡形态

由图16可以看出,空化器预置舵角为10°时,空泡轴线向上倾斜,航行体圆柱段仍包裹于空泡内。空化器预置舵角为20°时,空泡轴线上倾状态更加明显,航行体圆柱段出现沾湿(见图中红色区域),这是由于空化器预置舵角的存在,导致航行体周围的空化流场呈上下不对称分布,同时空化器法线与来流速度的夹角Δδ,相当于产生了附加攻角,空化器产生的升力会造成空泡轴线偏斜[15]。

图17给出了空化器预置舵角β为10°和20°,摆动频率f=1.5 Hz,航行体连续摆动过程中升力系数随摆动攻角的变化曲线。表3统计了摆动周期内的升力系数峰值与攻角对应关系。可以看出,首次负攻角摆动过程中,航行体受到正向的法向力冲击,量值与舵角呈正相关,随着预置舵角增大,升力冲击峰值增大,且出现的时间提前;随后的往复回摆运动过程中,升力系数为恒定值时,航行体处于空泡完全包裹状态,此阶段航行体受力与航行攻角的方向无关,表现为稳定的低头升力或抬头升力,这与袁绪龙等[16]获得的预置舵角能够产生稳定的抬头升力从而实现弹道转平的结论一致。在航行体往复回摆过程中,空化器预置舵角的存在,使得流体动力延迟特性在经历不同正/负攻角时呈现不同特性。正攻角运动阶段中,随着预置舵角增大,延迟效应减弱,甚至出现相位提前,升力峰值减小;负攻角运动阶段中,预置舵角导致延迟效应加剧,升力峰值增大。

表3 升力系数峰值与攻角对应关系

图17 升力系数随摆动攻角的变化曲线

为了进一步分析预置舵角对航行体流体动力延迟特性影响的形成原因,图18给出了图17中特征点1,2前后的空泡特征。由图18分析特征点1的形成原因。在t=1.030 s时刻,航行体只有空化器沾湿,受到由预置舵角形成的恒定负向法向力,并随着航行攻角的减小而减小。随着航行体下摆运动,发生空泡上飘现象,加之空泡的改变相对于航行体运动的滞后性,导致了航行体尾部背流面开始沾湿,航行体受到了较大的正向升力,见t=1.045 s时刻。随着航行体负向攻角的继续增大,空泡上飘现象逐渐被抑制,空泡轴线在流体作用下开始下偏,见t=1.060 s时刻。特征点2的形成原因分析:当航行体进入下一个运动周期时,空泡延迟效应产生的空泡轮廓弯曲导致圆柱段下方出现沾湿,且沾湿面积逐渐扩大,见t=1.701 s时刻。在t=1.705 s时,航行体圆柱段下侧沾湿面积达到最大后迅速降低造成了法向力曲线上的尖点,摆动到t=1.730 s时航行体被空泡完全包裹。

图18 不同时刻的空泡特征

3 结 论

本文采用数值仿真方法研究了航行体连续变攻角过程中非定常流体动力延迟特性,获得了速度、摆动频率和空化器预置舵角对航行体流体动力特性的影响规律与形成原因。得到的结论如下:

1) 航行体非定常超空泡航行状态下的流体动力具有延迟特性,相比于定常超空泡计算结果,自由攻角范围增大,同一航行攻角下的升力系数减小。

2) 非定常流体动力主要由航行体圆柱段上的非对称闭合产生,附体空泡的动态发展是形成升力非线性变化的原因。摆动频率的增加,会加快附体空泡与主空泡的融合,从而减弱自由空泡溃灭产生的冲击载荷,同时,对升力特性的延迟效应减弱。

3) 超空泡航行体在连续变攻角运动过程中,速度的增加使得附体空泡溃灭加剧,产生的冲击载荷最高可达主峰峰值的2倍;空化器预置舵角的存在导致航行体周围空化流场的不对称分布,在空泡延迟效应和空泡轴线偏移的共同作用下,航行体正负攻角摆动时具有不同的流体动力延迟特性。