孟旭飞 白 鹏 刘传振 李 盾 王 荣
(中国航天空气动力技术研究院,北京 100074)
随着各国军事水平及空间探索活动的进一步发展,各类飞行任务的复杂性及执行成本日益增加,各国对低成本、高效率运载器的需求更加旺盛.水平起降的吸气式动力高超声速飞行器作为一种大气层内高效运载手段,利用可重复使用特性可降低运载成本.目前,高超声速飞行器大多采用细长体、组合体、融合体和升力体气动布局形式,但面对水平起降所需的宽速域高升阻比的气动设计需求,还难以完全满足.
乘波体因其高超声速、高升阻比特性,目前已应用于多种高超声速飞行器布局选型中[1].乘波体迄今发展已有50 多年历史,1959 年Nonweiler[2]首先提出了楔形流场乘波体,1968 年Jones 等[3]使用锥形流场设计乘波体,有效提高了容积率,后续学者们研究使用了其他的流场,包括带攻角锥、椭圆锥流动[4-5]、一般三维流动[6]等,拓展了乘波体的设计空间.1990 年Sobiezky 等[7]提出了密切锥方法,大大提高了设计灵活性,可以得到具有特殊性质的乘波体外形[8-9].
但是,由于乘波体存在低速性能差、稳定性不好等缺陷,极大地限制了在工程中的应用.在提高低速性能方面,学者们做了一些努力.其中,为乘波体设定合理的平面形状是一种可行思路,李珺等[10]将双后掠平面形状投影到锥导激波流场生成外形,改善了低速性能;美国洛克希德马丁公司的Rodi[11]从密切锥方法出发,提出定后掠角密切锥乘波体的新概念[12];Zhao 等[13]根据Rodi 的思路设计了一种尖头后掠乘波体外形,不仅保持了高超声速性能,而且在涡升力作用下,低速性能也有所提高[14].刘传振等[15-16]提出过一种双后掠乘波体,发现外形在保证高超声速高性能的同时,低速性能较好.
而在改善乘波体的稳定性方面,相关研究还不多见.贾子安等[17]研究过乘波体的纵向稳定性,指出基于锥形流场的外形不利于提高纵向稳定性;Bykerk等[18]分析了乘波体飞行器低速状态的纵向稳定性,发现重心位置非常靠前时才能全程静稳定.同时,在乘波体横侧向稳定性[19]方面的工作更是少见.一般来说,上下反外翼是改变飞行器横侧向稳定性的有效思路[20-21],这在低速飞行器中比较常见,而He 等[22]基于锥形流场设计过外翼上反的乘波体外形,发现上反翼可以改变压心位置,有利于提高静稳定性.
因此,综合改善乘波体低速特性和稳定性的思路,为双后掠乘波体设计上/下反外翼,应该是一项值得研究的课题.本文根据密切锥方法提出了给定前缘线的乘波体设计,生成分别具有上反和下反外翼的双后掠乘波体,之前已有工作通过对理想化乘波体外形进行研究,验证双后掠乘波体保持了良好的高超声速特性[15],本文则关注低速状态,并重点研究使用上下反翼改善低速稳定性的可行性,为乘波体外形在宽速域高超声速飞行器设计中提供新的思路,进而提升总体性能.
密切锥乘波体方法由给定的激波出口形状拟合流场,进而生成乘波体外形.Sobieczky 等[7]提出的密切锥理论假定: (1) 每个密切面内流场的激波与指定的激波形状一致;(2)相邻密切面内的横向流动足够小.在密切锥方法中,有两条设计曲线,一条是激波出口型线(inlet capture curve,ICC),决定激波的出口形状;另一条为流线追踪起始线(flow capture tube,FCT),表示流线追踪的起始投影线.
在密切锥乘波体设计中,设计曲线ICC,FCT与俯视图上的平面形状轮廓线(planform contour line,PCL)存在几何关联.选择标准坐标系,以飞行器的展向y坐标为自变量,将ICC,FCT 和PCL 由3 个函数c(y),f(y)和p(y)表示,那么函数c,f和p之间的几何关系可用如下的微分方程组表达
式中βc为锥形流激波角,具体推导请参考文献[23].边界条件为3 条曲线的交点K,即f(yK)=c(yK)=p(yK),yK为半展长.
微分方程组(1)建立了密切锥乘波体的设计参数ICC,FCT 与平面形状PCL 之间的几何关系,我们命名为刘-白方程.观察此方程组,未知数有x,y两个,βc为已知量,根据微分方程理论只要已知函数f,c,p中的2 个,第3 个就可以求出.
存在3 种情况: (1) 已知f和c,求p;(2) 已知f和p,求c;(3) 已知p和c,求f.因为函数f和p分别表达前缘线在正视图和俯视图方向的投影,所以情况(2)可以看作给定三维前缘线设计乘波体外形的方法.
本节通过给定前缘线生成外翼上反和下反的乘波体.基准外形来自文献[15]中的弯头双后掠乘波体,长度8 m,翼展9 m,第一后掠角75°,第二后掠角50°.
提取双后掠乘波体的前缘线,如图1 所示的“basic”曲线,它在正视图方向为水平直线.保证弯头双后掠俯视图平面形状不变,生成具有上/下反外翼的前缘线,如图1 中所示,分别标记为dih-2,dih-1,anh-1,anh-2,4 条曲线.这4 条曲线与“basic”有相同的俯视图投影,而在正视图方向,翼梢分别抬高0.8 m,0.4 m 和降低0.4 m,0.8 m.
图1 基于双后掠乘波体提取的前缘线Fig.1 Leading edges from double swept waverider
由这些前缘线生成乘波体外形,所有曲线均由离散点形式表达,每条曲线包含201 个离散点.设定锥形流激波角βc=15°,设计状态H=30 km,Ma=5.将前缘线沿俯视图和正视图方向投影得到p和f,使用数值方法求解方程(1)得到函数c.图2 左侧给出根据这5 条前缘线求解得到的ICC 曲线,右侧则是将ICC 和FCT 输入密切锥方法生成的乘波体外形.可以看到,dih-2 和dih-1 外形为外翼上反外形,dih-2 的上反幅度更大;anh-1 和anh-2 为外翼下反的外形.这5 个外形具有相同的俯视图平面投影形状,并且容积率相同,均为0.141 8.
图2 求解的ICC 曲线和对应外形Fig.2 Solved ICCs and generated waveriders
图2 求解的ICC 曲线和对应外形(续)Fig.2 Solved ICCs and generated waveriders (continued)
使用数值模拟手段评估外形的气动性能.数值方法采用有限体积方法求解三维可压缩N-S 方程.无黏通量采用Roe 格式[24]计算,有权重格林-高斯公式重构方法获得空间二阶精度,梯度限制器选取改进的Barth 限制器[25],以消除计算中间断附近的数值过冲和振荡,同时应用基于压力辅助限制器的局部熵修正方法避免非物理解[26].黏性通量采用二阶中心格式计算.湍流模型采用在工程上广泛应用的Menter SSTk-ω两方程模型[27].时间方向采用二阶精度双时间步方法,LU-SGS 隐式推进求解[28].计算平台为我院自主研发的计算流体力学(computational fluid dynamics,CFD)数值模拟软件平台—GiAT,已经过大量理论和工程验证[29-32].
为了保证黏性流动的精确模拟,采用分区结构化网格.黏性层法向网格尽量保证垂直壁面,预估第一层边界层高度保证y+量级为1,流向参数梯度较大的位置保证足够的流向网格密度,计算网格如图3 所示.
图3 计算网格Fig.3 Computational mesh
以Basic 外形为例分析网格收敛性.分别生成3 套不同尺度的网格: 稀疏网格Coarse,网格数422 万左右;中等网格Medium,网格数818 万左右;密网格Refined,网格数1610 万左右.
在H=0 km,Ma=0.4 状态评估气动性能,首先考察升阻力数据.表1 给出在攻角α=10°,侧滑角β=0°,力矩参考点取头部顶点时升阻特性计算结果,可以看到Medium 网格与Refined 网格之间差异较小,均在0.3%以内,而Coarse 网格则与密网格差异较大,最大差异为升力系数1.67%.
表1 不同网格的升阻特性结果Table 1 Lift-drag results of different grids
表2 给出力矩参考点取头部顶点,α=10°,β=3°时滚转力矩、偏航力矩和俯仰力矩的结果,同样可以看到Medium 网格与Refined 网格之间差异较小,均在0.8%以内,Coarse 网格与Refined 网格最大差异为偏航力矩3.36%.
表2 不同网格的力矩结果Table 2 Moments results of different grids
对于乘波体横侧向稳定性导数及侧向力关于侧滑角的导数对网格敏感性研究结果在表3 中给出.以密网格作为基准可知,Medium 网格与Refined 网格之间差异更小,变化范围均在0.8%以内,Coarse网格与Refined 网格差异较大,最大差异达到6.17%.
表3 不同网格稳定性导数结果Table 3 Lateral-directional derivatives results of different grids
考虑到Medium 网格相比Refined 在气动力、力矩和稳定性导数方面的误差较小,而Coarse 网格误差较大,本文认为Medium 网格,即网格量818 万左右,是可信和足够的.
本节研究上/下反翼对双后掠乘波体升阻特性的影响.
计算图2 中5 个外形的气动力,计算状态为H=0 km,Ma=0.4,α=0°~ 20°,β=0°.这5 个外形具有相同的平面投影形状和容积,参考面积为29.598 m2,参考长度为8 m.
图4 给出5 个乘波体升阻力特性随攻角的变化.可以看到,带上/下反翼乘波体与basic 外形的升阻特性变化规律一致,最大升阻比均在α=4°附近取得.相比于basic 外形,机翼上反使得乘波体升力和阻力均有所减小,机翼下反使得升力和阻力增大,综合导致在α=4°时anh-2,dih-2 升阻比较basic 外形分别高出0.77% 和3.48%,anh-1,dih-1 升阻比较basic 外形分别降低0.76%和2.74%.通过与高超声速状态对比可知,机翼上下反对乘波体升阻特性影响在低速阶段更为明显,上下反程度较大时对升阻比有一定的提升,但幅度较小.
图4 升阻特性随攻角的变化(有底阻)Fig.4 Lift-drag characteristics variation with angle-of-attack (with bottom drag)
低速状态下底阻对飞行器气动性能影响较大,实际设计中需要对后体进行修形减阻,同时此类宽速域飞行器一般将发动机尾喷口置于飞行器底部,底阻一般小于本文中所采用外形,因此在实际设计评估中需要对底阻扣除或进行相应处理.
图2 中的5 个外形均由密切锥乘波体设计方法生成,具有相同的平面投影形状,2.2 节的气动力分析也说明这几个外形具有较为一致的升阻特性变化规律,同时不同外形间升阻特性存在一定差异.本节通过流场结构分析,进一步说明上下反翼对乘波体低速升阻特性的影响规律.
低速条件下,乘波体背风面旋涡会对气动性能产生影响.图5 给出图2 中不同外形在最大升阻比点(α=4°,β=0°)时,上表面压力分布及涡结构,可以看到,这5 个外形背风面均具有明显对称涡结构,造成机翼两侧大范围的低压区,使得乘波体低速升力特性得到改善,上表面压力分布近似,不同外形间升阻特性差异较小.对于basic 外形,下表面和上表面阻力占比分别为23.43%和26.94%,升力占比分别为25.35%和75.40%,且不同外形上下表面升阻力占比相近.图6 给出了不同外形间下表面压力分布,可以看到不同外形间差异主要体现在第一后掠角和第二后掠角交接处沿前缘向后延伸的高压区域范围,机翼下反使得高压区范围增大,下表面提供的升力和阻力均增大,反之,机翼上反使得下表面提供的升力和阻力均减小,综合作用下对升阻特性影响结论如2.2 节所述.
图5 上表面压力分布及涡结构Fig.5 Pressure distributions of upper surface and vortex structure
图6 下表面压力分布Fig.6 Pressure distributions of lower surface
图6 下表面压力分布(续)Fig.6 Pressure distributions of lower surface (continued)
本节研究外翼上/下反对纵向稳定性的影响.一般飞行器设计中将重心布置于2/3 全长处,图7 给出了力矩参考点取66.7%全长处时俯仰力矩随攻角的变化.可以看到,不同外形的俯仰力矩变化趋势比较类似,同时也存在一定的变化差异,且随攻角增大而增大.在α=20°时,相比于basic 外形,dih-2 的俯仰力矩系数增加11.92%,dih-1 增加8.82%,anh-1 减小6.96%,anh-2 减小13.08%.
图7 俯仰力矩系数随攻角变化Fig.7 Pitching moment variations with angle-of-attack
图8 给出了不同外形气动焦点(aerodynamic center,A.C)随攻角的变化.纵向稳定性可以根据气动焦点与重心的位置判定,当重心在气动焦点之前时,飞行器是纵向稳定的.由图8 可知,当α> 4°时各外形气动焦点都位于66.7%之前,所有外形均为纵向静不稳定.机翼上反时,气动焦点前移,纵向稳定性降低;机翼下反时,气动焦点后移,纵向稳定性提高,该影响规律与高超声速条件下上下反翼对乘波体纵向稳定性影响规律一致,因此,机翼下反有助于提升乘波体全速域的纵向静稳定性.
图8 气动焦点随攻角变化情况Fig.8 Pressure centers variations with angle-of-attack
当然气动焦点位置变化量较小,说明纵向稳定性主要跟平面形状有关,上下反外翼影响有限.
目前常用的评估横侧向静稳定性的判据有[33]:Clβ,Cnβ,Cnβ,dyn和LCDP,其中Clβ和Cnβ表示横向与侧向的单通道静稳定性,Cnβ,dyn为偏航动态稳定判据,表征飞行器开环无控状态横侧向稳定特性.LCDP 为横向控制偏离判据,用于预测闭环有控状态横侧向控制稳定性,这种形式的横侧向控制偏离判据与控制策略相关,控制策略不同判据表达式也不同,且求解过程涉及气动导数过多,计算复杂.
本文采用形式较为简单的Clβ,Cnβ和Cnβ,dyn作为稳定性判据评估上/下反翼对乘波体高超声速横侧向稳定性的影响.需要指出的是,本文采用常用的机体坐标系,坐标系定义为:X轴从质心指向飞行器头部,Y轴从质心指向上,Z轴由右手定则确定.在这种坐标系下,Clβ< 0 表示横向静稳定,Clβ> 0 表示横向静不稳定;Cnβ< 0 表示侧向静稳定,Cnβ> 0 表示侧向静不稳定;Cnβ,dyn< 0 表示开环状态偏航通道动态稳定,Cnβ,dyn> 0 表示开环状态偏航通道动态不稳定.
为了验证Cl,Cn在计算Clβ,Cnβ的侧滑角范围内与β为线性相关,首先给出了β=0°~ 8°范围内Cl,Cn随β的变化曲线,如图9.可以看到Cl,Cn与β线性相关性较好,基于此Clβ,Cnβ只需选取两个固定侧滑角状态即可得到.
图9 线性验证Fig.9 Linearity verification
先分析横向静稳定性.图10 给出了不同外形Clβ随攻角的变化,可以看到不同外形的Clβ均随攻角增加呈先减小后增大的趋势.同一攻角时,对比basic 外形,机翼上反使得Clβ减小,且Clβ随机翼上反程度增大而减小,机翼下反使得Clβ增大,这说明机翼上反有助于提升乘波体的横向静稳定性,而机翼下反会使横向静稳定性下降.
图10 Clβ 随攻角的变化Fig.10 Rolling-moment derivative variations with angle-of-attack
侧向静稳定性方面,图11 给出Cnβ随攻角的变化,可以看到,5 种外形乘波体Cnβ均小于0,即均为侧向静稳定的.对比basic 外形,机翼上反的外形Cnβ负值更大,即侧向静稳定性较basic 有所提高,机翼下反使得侧向静稳定性降低.dih-2 外形Cnβ负值更大,说明机翼上反程度更大,乘波体侧向静稳定性更好.
图11 Cnβ 随攻角的变化Fig.11 Yawing-moment derivative variations with angle-of-attack
图12 给出了α=10°,β=10°,上表面压力分布及流场涡结构示意.可以看到,在侧滑角不为0°时,上表面两侧低压区存在明显差异,不同外形间两侧低压区压力分布位置相近.通过不同切面流场涡量图也可以看到,不同上下反翼外形流场涡结构一致,由头部第一后掠区向后延伸,且左侧涡明显强于右侧,对应上表面两侧低压区形态.不同外形间左侧低压区范围及压力分布较为一致,差异主要体现在右侧低压区范围及翼梢处的压力分布.相比于无反外形,机翼上反使得右侧低压区明显减弱,两侧低压区差异增大,翼梢处高压区范围增大,综合导致滚转力矩增大,横向稳定性提升;机翼下反使得右侧低压区明显增强,两侧低压区差异减小,翼梢处高压区范围差异较小,综合导致滚转力矩减小,横向稳定性降低;图13 给出了不同外形α=10°,β=10°时不同外形下表面压力分布,可以看到各外形间下表面压力分布近似,机翼下反时右侧翼梢处出现小范围高压区,使得两侧侧压力差异减小,滚转力矩减小,横向稳定性下降.
图12 涡结构示意Fig.12 Vortex structure
图13 下表面压力分布Fig.13 Pressure distributions of the lower surface
下面考察偏航动态稳定性判据Cnβ,dyn,具体表达式[34]为
其中Ixx和Iyy分别为飞行器横向转动惯量和航向转动惯量,本文中以basic 外形为基准外形,假设该外形为简单的铝制壳体结构,调整内部空间使质心位置位于2/3 全长处,此时测得相对质心转动惯量Ixx/Iyy为2.938.
由公式(2)可知,随攻角增大,Clβ对偏航动态稳定性贡献度增大,由于横向静稳定性裕度随攻角增大而增大,此时飞行器偏航动态稳定性得到提升.图14给出了Cnβ,dyn随攻角的变化情况,可以看到,在计算状态下,所有外形均为偏航动态稳定状态.对比basic外形,机翼上反使乘波体偏航动态稳定性明显提升,且提升程度随机翼上反程度增大而增大;机翼下反使乘波体偏航动态稳定性降低,且降低程度随机翼下反程度增大而增大.
图14 Cnβ,dyn 随攻角的变化Fig.14 Lateral-directional dynamic stability variations with angle-of-attack
本文基于密切锥方法发展了给定前缘线的乘波体设计,通过该方法生成理想的尖化前缘乘波体外形,研究上下反翼对双后掠乘波体低速特性的影响,发现:
(1) 低速状态下机翼上下反对乘波体升阻特性影响较高超声速条件下更加明显,但总体影响有限,仅在机翼上下反程度较大时升阻比有一定提升;
(2) 纵向稳定性方面,低速状态下各外形为静不稳定,这与高超声速状态不同,设计时需注意.机翼上反,乘波体气动焦点前移,纵向稳定性降低;机翼下反,气动焦点后移,纵向稳定性提高;
(3) 横侧向稳定性方面,机翼下反使得乘波体下表面背风区翼梢处产生小范围的高压区,减小了下表面两侧压力差异,导致滚转力矩减小,横向静稳定性下降,机翼上反则使上表面两侧压差增大,横向静稳定性提升;同样的,机翼上下反所引起的两侧压力差异变化导致偏航力矩大小发生改变,机翼上反则使上表面两侧压差增大,侧向静稳定性提升;
(4) 偏航动态稳定性方面,在横侧向静稳定性综合作用下,机翼上反使乘波体偏航动态稳定性有明显提升;下反则会降低.