基于负压梯度翼型栅格翼设计与水中数值分析

2023-04-06朱滢铧邹志辉蒋运华

弹道学报 2023年1期

朱滢铧,邹志辉,李洋,沈乐,蒋运华

(中山大学海洋工程与技术学院,广东珠海 519000)

水下发射的航行体出水过程涉及复杂的多相流动,常常出现俯仰角速度过大或弹体结构破坏等问题,从而导致出水失败。栅格翼应用于水下航行体的设计,有望提高其出水稳定性。由于负压梯度外形栅格翼具备较好的流动分离特性,能有效推迟空化发生,有望应用于水下栅格翼的设计中。

NESTERUK[1-2]研究表明较短的正压梯度区的外形可以改善减阻性能。BURAGA等[3]比较了超空泡流和不分离流两种流态下的轴对称细长体的阻力,结果表明空化数小于0.001时,超空泡流几乎没有优势。WU等[4]计算了不同构型栅格翼在不同空化数和攻角下的水动力特性,结果表明,当流动条件相同时栅格翼的升力并不会随着叶片数量的增加而增加。苏太昌等[5]研究了不同剖面形状的栅格及不同空化数对栅格翼升阻力系数的影响,结果表明,当空化数相同时尖锐的栅格剖面形状可以明显提高升力。夏艳艳等[6]对不同构型栅格翼的水动力特性进行研究,结果表明增加栅格翼叶片数量会同时增加升力和阻力,发生空化后的栅格翼升力会急剧下降。李磊[7]设计了一种具有较长的逆压区间和较小的逆压梯度的轴对称细长体模型来实现无分离流动。郭小飞[8]针对通过使用不分离外形来降低压差阻力的减阻方式,做了数值模拟的侧面验证,结果表明这一水下减阻方案在一定速度范围内是可行的。安伟光等[9]通过数值模拟对负压梯度轴对称外形的分离特性进行了探究,验证了势流理论与NESTERUK[10]风洞试验的压力分布,并对外形进行空化计算,结果表明负压梯度轴对称外形发生了空化现象,并且形成了超空泡。

栅格翼应用于水下航行体时会导致阻力大幅度增加并引起空化现象,尤其是在攻角较大时,空化的发生会使升力急剧下降,从而导致控制效率降低。因此,设计水下易空化下具有优良升阻特性的栅格翼是其应用的关键。本文依据势流理论,设计了具有较长的负压梯度区的流动外形,并以此为基础设计了几种框架型的栅格翼。通过数值模拟验证了不分离特性并最终选出最适合用于航行体姿态稳定的栅格翼布局,可为负压梯度外形的理论研究和新型栅格翼的设计工作提供一定的参考。

1 计算模型和工况及网格划分

基于负压梯度外形曲线方程[9],设计了一种负压梯度(UA)外形,其弦长为L,厚度为R,最大厚度点x*位于三分之一弦长处(L/3)。其外形曲线与压力分布如图1所示。图中,Cp为翼型表面的压力系数。

图1 UA不分离流动外形与其压力分布Fig.1 The negative pressure gradient foil and its pressure distribution

图2为基于UA外形的3种不同布局的栅格翼。图2(a)为三维情况下的一种栅格翼布局,图2(b)为所设计的3种布局。每种布局下的UA截面尺寸一致。3种布局的d/L(d为叶片间距)分别为1.71,1.11和0.85。

图2 3种栅格翼布局Fig.2 Three kinds of structural layouts of grid fins

图3(a)和图3(b)分别为UA和新型栅格翼的计算域与网格划分。对于UA的计算工况：以翼型弦长为特征长度的弗汝德数Fr=6.57,空化数σ=1,攻角α=0°,7°。采用结构化网格对计算域进行划分,在水翼来流区、尾迹区和翼型周围进行加密。对于新型栅格翼的工况：Fr=65.36;σ=0.2;α=-8°,-6°,-4°,-2°,0°,2°,4°,6°,8°。计算域尺寸与UA保持一致,依旧采用结构化网格划分。

图3 UA和新型栅格翼的计算域与网格划分示意图Fig.3 The computing domains and grid divisions

2 数值方法与有效性验证

本文多相流模型采用Mixture模型,湍流模型采用DES模拟,在该模拟中近壁面的RANS采用SSTk-ω模型,空化模型采用Zwart-Gerber-Belamri模型。为验证网格无关性,选取同尺寸的NACA0015二维翼型,计算工况：Fr=6.57,α=7°,σ=1。本文采用了5套网格进行计算与比较。由于计算中均使用了壁面函数,其中3套网格的第一层Y+都保证处于30～60之间,网格量分别为216 130,512 980和803 090,分别记为网格Ⅰ、网格Ⅱ与网格Ⅲ。此外,考虑到本文算例对于翼型表面水动力参数的计算精度要求,增加网格Ⅳ(网格量为512 980,Y+=1)与网格Ⅴ(网格量为803 090,Y+=1),网格Ⅳ中翼型的网格划分如图4所示。

图4 NACA0015翼型的网格局部放大Fig.4 The grid of the NACA0015 airfoil

将得到的升力系数CL和阻力系数CD与文献[11]的实验数据相比较,对比结果如表1所示。可以看出,网格Ⅳ的升力系数精度与阻力系数精度都要明显高于网格Ⅰ、网格Ⅱ和网格Ⅲ。网格Ⅳ的升力系数精度与网格Ⅴ的接近,虽然网格Ⅳ的阻力系数精度不如网格Ⅴ,但考虑到计算时间成本,最终选取网格Ⅳ进行后续的计算。

表1 不同网格数下的升阻系数验证Table 1 Lift and drag coefficients under different grid sizes

图5为网格Ⅳ压力系数的计算结果与文献[12-13]中实验数据的对比图(Fr=6.57,α=7°,σ=1)。可以看到,计算结果与文献[12]的实验数据基本吻合,与文献[13]的实验数据有较高的一致性。

为进一步验证不分离外形的流动特性,对UA外形进行了数值模拟。图6给出了Fr=6.57,α=0°,σ=1时UA外形表面的壁面剪应力(τ)分布曲线。可以看出,UA外形的壁面剪应力在整个表面均大于0,所以并未出现流动分离,验证了设计理论的有效性。

图6 UA外形表面的壁面剪应力分布曲线Fig.6 The wall shear stress distribution of UA profile

3 计算结果与分析

3.1 负压梯度翼型的水动力特性

为探究UA外形的水动力特性并验证其不分离优势,将其计算结果与同种工况下NACA0015翼型的水动力特性进行对比。图7给出了Fr=6.57,α=7°,σ=1下UA与NACA0015翼型的升阻系数对比;图8(a)和图8(b)分别给出了该工况下NACA0015翼型和UA翼型在一个周期内空泡脱落的水的体积分数图。

图7 UA与NACA0015翼型升阻系数对比Fig.7 The lift and drag coefficients of UA profile and NACA0015 airfoil

从图8可以看出,在考虑空化流动后,二者均出现明显的空泡脱落,脱落频率和体积分数均表现出周期性变化。随着空泡脱落,二者的升力系数、阻力系数无明显差异,均表现出强烈的周期性。

图8 NACA0015与UA翼型一个周期内空泡脱落的水的体积分数图Fig.8 The shedding of cavities in one period of NACA0015 airfoil and UA profile

图9为Fr=6.57,α=7°,σ=1时UA与NACA0015翼型的压力系数分布曲线。临界空化数定义为

图9 UA与NACA0015翼型的压力分布Fig.9 The pressure coefficent distribution of UA profile and NACA0015 airfoil

式中：p∞为远场压力;pv为饱和蒸汽压力;ρ为流体密度;vcr为临界速度,本文中为外部主流速度。理论上临界空化数越小,也即最小压力系数的绝对值越小,越不易发生空化。

可以看出,NACA0015翼型压力面上压力系数分布为顺压梯度—逆压梯度;而UA外形则为顺压梯度—逆压梯度—较长一段区域的顺压梯度。这种压力分布会比传统的顺压梯度后只有逆压梯度的压力分布要更有优势。对于吸力面上的压力分布,UA外形吸力面上所有点的压力系数基本都比NACA0015要更靠近Cp=0。而从整体上看,UA外形上|Cp,min|比NACA0015的要小,也即其对应的临界空化数要比NACA0015的小。因此,UA外形的压力分布在水下应用时预计会比NACA0015翼型更有优势。

3.2 负压梯度翼型栅格翼的水动力特性

3.2.1 叶片间距对负压梯度翼型栅格翼水动力的影响

图10给出了3种布局下栅格翼升力系数和阻力系数随攻角的变化关系。对于升力系数关系,可以看出,3种布局下的|CL|随|α|变化趋势基本相同,|CL|均随|α|增大而增大。此外,叶片数量越多(叶片间距越小),|CL|随|α|的变化速率越快。因而导致,随|α|增大不同布局下的升力系数差值|ΔCL|明显增大,正负攻角下差值基本相等,参见图11给出的各攻角下d/L=0.85和d/L=1.71两种布局的升阻系数差值。

图10 3种布局下栅格翼的升、阻力系数随攻角的变化情况Fig.10 The lift and drag coefficients of grid fins with three kinds of structural layouts under different attack angles

由图10(b)可以看出,3种布局下的|CD|随|α|变化趋势也基本相同,|CD|均随|α|增大而增大,CD-α整体曲线变化基本关于α=0°曲线对称。对于任意攻角,CD均随叶片数量增多而增大,具体原因于3.2.2节中进行探讨。此外,随|α|增大,不同布局下的|ΔCD|也增大,但相对升力系数来说差值增大得不明显,且正负攻角下差值基本相等,具体如图11所示。

图11 各攻角下d/L=0.85和d/L=1.71的升、阻系数差值Fig.11 The lift and drag coefficients differences of grid fin with d/L=0.85 and grid fin with d/L=1.71 under different attack angles

当|α|较大时,|ΔCL|与|ΔCD|接近,但随着|α|变小,|ΔCL|急剧减小,|ΔCD|值相对|ΔCL|来说基本不变,如图11所示。因此,当攻角较小时,随叶片数量增加,CL值趋于一个常值,而CD值增长较为剧烈。此时,为保证具有良好的升阻特性,栅格翼的叶片数量在保证升力条件下应尽量少。当攻角较大时,CL值和CD值均随叶片数量增加而增大。此时无法直接判断栅格翼升阻特性的变化情况,这就需要研究各个布局下的升阻比特性。

图12 3种布局下栅格翼的升阻比随攻角的变化曲线Fig.12 The lift-drag characteristics of grid fins with three kinds of structural layouts under different attack angles

3.2.2 负压梯度翼型栅格翼的流动特性

由3.2.1可知,对于任意攻角,阻力系数均随叶片数量增多而增大。这是由于叶片间距越小,各叶片的压力分布相互干扰越剧烈,最终使得压差阻力增加。图13给出了0°、4°和8°攻角时3种布局下的栅格翼压力分布云图。可以看出,无论何种布局以何种攻角运动,从整个压力分布看,高压区均基本聚集在翼型前缘处,而翼型后部压力均小于前部压力,从而形成压差阻力。但是,翼型前部高压区的分布情况明显不同,使得不同布局下的压差阻力不同。

图13 当攻角为0°、4°和8°攻角时3种布局下的栅格翼压力分布云图(Fr=65.36,σ=0.20)Fig.13 The pressure distribution of grid fins with three kinds of structural layouts under the attack of angles of 0°,4° and 8°

对于d/L=0.85布局(图13(a)～图13(c)),此时翼型前缘高压区的压力均产生了强烈的干扰,这就导致栅格翼整体上前部的高压区分布十分明显,使得压差阻力明显增大。而d/L=1.11布局(图13(d)～图13(f))下的压力干扰相比d/L=0.85布局的情况有所减弱,这种布局下的压差阻力比d/L=0.85布局的小,但数值依旧较大。对于d/L=1.71布局(图13(g)～图13(i)),此时栅格翼前缘的压力干扰均很小,整体上翼型前缘的高压区分布不显著,前后压力差不明显,压差阻力也相对较小。所以,任意攻角下阻力系数随叶片数量增多而增大。

栅格翼前端的压力干扰会显著影响其水动力特性,图14给出了0°、4°和8°攻角时3种布局下的水的体积分数分布图。

对于同一布局、不同攻角的情况,每种布局下攻角由0°到4°再到8°变化时,各叶片空泡形态会发生变化。0°攻角下(图14(a)、图14(d)和图14(g)),各叶片空泡形态关于栅格翼轴线对称,最上方与最下方叶片的空泡长度最小,中轴线处叶片的空泡长度最大,但厚度无明显差别。4°攻角下(图14(b)、图14(e)和图14(h)),由上至下空泡长度与厚度均减小。这是由于,当栅格翼以一定攻角迎接来流时,各叶片迎流面是叶片下表面,此时下表面是压力面。而当前叶片下端的叶片则处于高压区流动,这使得两叶片间流区中的流线弯曲,本质上减少了当前叶片下端叶片的当地攻角。从而导致空泡形态自上而下逐渐变短、变薄,此外,空泡同样受重力影响,这可以理解为多重因素造成的结果[14]。而8°攻角下(图14(c)、图14(f)和图14(i)),各叶片空泡形态分布与4°情况基本一致,区别在于由上至下空泡长度与厚度的减小速率更大。

图14 当攻角为0°、4°和8°攻角时3种布局下的水的体积分数分布图Fig.14 The volume fraction of water of grid fins with three kinds of structural layouts under the attack of angles of 0°,4° and 8°

对于同一攻角、不同布局的情况,当栅格翼处于一定攻角时,随着叶片间距的增加,其空泡长度逐渐变小。这是因为,当叶片间距较大时(图13(g)～图13(i)),各叶片之间的压力相互干扰程度较其他两种布局的小,接近于单个叶片的流动特性。而当叶片间距减小时,整个栅格翼周围的压力均受到强烈的干扰,低压片区也相对较大,所以此时空泡形态更长。