0 引言
巡飞弹是一种低成本的精确制导弹药,具有长时间续航能力,并能够快速攻击陆地或海上非视线下的目标[1]。由于战备需求,巡飞弹常处于低雷诺数环境[2]。然而,目前有关翼型的研究大多基于传统飞行器翼型设计,且主要关注高雷诺数(Re=107量级)下的翼型优化。雷诺数从低到高所反映的气动现象完全不同,若直接将高雷诺数翼型套用进低雷诺数(Re=104～105)巡飞弹将导致气动性能急剧恶化。基于这一背景,有必要对巡飞弹低雷诺数翼型进行专门的研究。
目前,针对低雷诺数翼型的气动优化设计相较于传统飞行器翼型并不是十分丰富,主要的研究手段为风洞试验和数值模拟。国内外尽管对翼型低雷诺数试验开展了大量的研究工作[3-4],但由于低雷诺数条件下流动稳定性差,风洞实验易受外界各种干扰,导致无法精确控制风洞流场环境,从而不同风洞实验获得的实验数据差异较大。因此常采用数值模拟对翼型作先一步优化研究,避免风洞实验的反复调试。
随着计算机性能的提高和数值方法的发展,对翼型优化设计问题的研究趋于采用翼型参数化结合优化算法模拟得出最优翼型的方向发展。廖炎平等[5]分别对对称翼型、弯度翼型和超临界翼型的几何外形进行最小二乘拟合,对比了型函数线性扰动法、特征参数描述法、正交基函数法和CST方法(class function/shape function transformation)对翼型参数的拟合精度,结果表明CST方法具有更好的拟合精度。自此,国内专家学者在翼型参数化时逐渐采用CST参数化作为最优方案。在优化算法方面,选择序列二次规划法[6],与其他算法相比,序列二次规划法求解非线性规划问题更加高效,既能保持整体收敛性又能完成局部超一次收敛性,更加适合实际工程需求[7]。
为提升巡飞弹在巡航阶段的气动性能,本文中采用CST参数化方法搭建翼型模型,以提高升阻比作为优化目标,根据数值模拟的气动力数据,结合序列二次规划法进行寻优,为低雷诺数巡飞弹翼型提供理论依据和设计参考。
1 CST翼型参数化模型建立
1.1 CST参数化方法
CST参数化方法(class function/shape function transformation)是美国波音公司的Kulfan等提出的一种基于型函数/类函数变换的参数化方法[8-9]。翼型的上下表面用如下的公式来表达:
上表面
(1)
下表面
(2)
式中:ψ=x/c为翼型无因次x坐标,ζ=z/c为翼型无因次y坐标;ΔξUpper=zuTE/C为翼型上表面前缘厚度比,ΔξLower=zlTE/C为翼型下表面后缘厚度比。
类函数
定义如下:
(3)
型函数S(ψ)定义如下:
(4)
(5)
(6)
式中:对于一般的圆头部、尖后缘翼型,N1取0.5,N2取1.0;Aui和Ali为翼型上下表面部件形函数的系数;Si(ψ)为 Bernstein多项式;n表示翼型上下表面Bernstein多项式的阶次。
因此,用CST参数化建模方法表示上下表面坐标方程表示为:
(ζ)Upper=
(ζ)lower=
1.2 不同Bernstein多项式阶数下精度对比
研究表明,Bernstein多项式阶数对CST参数化拟合精度影响较大[10-11]。Ceze对36阶Bernstein多项式定义的翼型参数化进行研究[12],发现过高的阶数将使参数化过程严重病态化,且不同翼型形状差异较大,针对亚音速翼型前缘半径较小,最大厚度位置靠后,低雷诺数翼型前缘厚度更厚,弯曲程度更为明显,因此不同翼型所需阶数并不相同,选取适当的阶数是合理且必要的。
结合国内外已知低雷诺数机翼研究,选取E387翼型作为巡飞弹基础翼型,弦长为100 mm,根据文献[13],为了保证上下表面在前缘处的前缘半径相同,令Au0=-Al0,为研究不同Bernstein多项式阶数对E387翼型参数化精度,取n=1～12进行研究。
图1为分别使用n=1～12时CST方法拟合翼型函数曲线点和实际翼型点的坐标残差的平方和。从图1可以看出,对E387翼型进行参数化表示时,Bernstein多项式阶数越高,误差减小,n>6时精度基本达到仿真要求,8阶以后,误差趋于稳定,可以达到很高的精度。图2为n=8时的CST参数化翼型与E387翼型对比,由图2可看出2种翼型基本重合,前缘和后缘翼型捕捉清晰,翼型弯度能很好满足,考虑到后续优化因变量不宜过多,因此本次优化仿真采用n=8作为Bernstein多项式阶数。
图1 不同阶数下残差平方和
Fig.1 Sum of squares of residuals for different orders
图2 n=8时,CST拟合结果和基础翼型对比
Fig.2 n=8, comparison of CST fitting results and basic airfoil
2 数值计算
2.1 网格划分
采用ICEM软件对设计翼型进行网格划分,采用C型拓扑结构,整个计算域可视为半圆形网格和矩形网格结合,翼型前端与尾部点到计算域边界的距离取E387翼型弦长的20倍,弦长为100 mm,第一层网格数取0.01 mm,翼型表面节点布置为330个,网格数量为48 960,流场网格划分整体轮廓图和局部加密图如图3所示。
图3 翼型网格划分图
Fig.3 Airfoil meshing diagram
2.2 控制方程及湍流模型
湍流流动的应用模型包含单方程模型、标准k-ω模型、修正的k-ω模型、大涡模拟等[14]。其中SST k-ω模型在近壁面采用k-ω模型,在边界层外采用k-ε模型,将两者优势相互结合,应用最为广泛。
利用Fluent压力-速度耦合求解器求解雷诺平均的不可压缩N-S方程,采用couple算法,梯度差值方法选用最小二乘法,压力差值方法选用二阶迎风格式,计算步长为200。
2.3 边界条件
翼型表面设置为无滑移翼面,进口处设置为速度进口,入流攻角为4°,速度为30 m/s,出口处设置为压力出口,计算得出升力系数和阻力系数。
2.4 验证案例
为验证本文中所采用的数值计算方法的可靠性,选择对E387翼型进行数值模拟,采用3.1节相同的网格划分方法,湍流模型分别采用S-A、SST k-ω、Transition SST,攻角分别设为-1.10、-0.09、1.01、1.95、3.00、4.03、5.02、6.04、7.08、8.12、9.07、10.08、11.22及12.14共13个攻角,并将仿真数据与伊利诺伊州立大学(UIUC)低速风洞在雷诺数Re=2×105条件下获得的E387翼型实验数据进行对比[15],对比结果如图4所示。
图4 3种湍流模型仿真数值与实验数据对比
Fig.4 Comparison of simulation values and experimental data of three turbulence models
从图4可以看出,SST k-ω、Transition SST相较于S-A湍流模型更加精确,与实验数据更为贴合。实际在仿真计算过程中,当攻角为4.03时,升力系数实验值为0.771,SST k-ω模型下获得的升力系数为0.788 66,偏差为2.29%,Transition SST获得的升力系数为0.788 79,偏差为2.31%;阻力系数实验值为0.013 5,SST k-ω模型下获得的阻力系数为0.012 81,偏差为5.11%,Transition SST获得的阻力系数为0.012 78,偏差为5.33%,相比而言SST k-ω模型略优于Transition SST模型,验证了本文中采用的数值计算方法的可靠性。
3 采用序列二次规划算法的气动优化设计
求解约束非线性优化问题时,序列二次规划算法是最有效的方法之一,具有收敛性好、计算效率高、边界搜索能力强等优点。
以E387翼型为基础翼型,要求设计一个低雷诺数翼型,使其在马赫数M=0.088,攻角α=4°,巡飞高度h=150 m,空气密度ρ=1.208 1 kg/m3,雷诺数Re=2×105时具有最大升阻比,同时翼型面积基本保持不变。采用Bernstein多项式阶数为8,则Aui和Ali(设计变量)共有18个。
每组新得出的18个变量由Isight循环到Matlab模块得出新的翼型参数坐标,通过编写的rpl脚本文件,在ICEM中自动划分网格。导出网格文件后,根据事先编好的jou文件,Fluent自动读取Mesh文件,完成对翼型的流场计算,得出升力系数和阻力系数。计算器模块完成升阻比的计算,优化后自动清除由Fluent导出的升力系数和阻力系数文件,以免下次循环时重复分析上次得出的气动结果。用序列二次规划算法进行寻优完成优化,流程如图5所示。
图5 采用序列二次规划算法优化设计流程图
Fig.5 Optimized design flow chart using SQP algorithm
在E387翼型优化设计阶段,给定翼型优化目标函数,其表达式为:
式中,以巡飞阶段4°攻角的升阻比最大值为目标,同时设置翼型最大厚度和最小厚度,来避免机翼在优化过程中可能发生的变形病态化。
根据上述参数设定,设置最大迭代次数为60,得出优化后的翼型结果。分析对比后数值如表1所示,可知翼型优化后的气动性能得到了提升,升阻比提高了1.917 3,提高百分比为4.173%,升力系数提高了0.060 19,提高百分比为11.127%,但同时阻力系数也升高了6.678%。
表1 优化前后气动力系数值
Table 1 Aerodynamic coefficients before and after optimization
模型ClCdCl/CdΔCl/Cd基础翼型0.540 950.011 7745.943 0优化翼型0.601 140.012 5647.860 44.173
4 气动特性分析
4.1 优化翼型气动分析
图6给出优化后翼型与基础翼型的几何形状相比,前缘厚度明显增加,后缘变化幅度较小,可以满足翼型在实际使用中的强度要求,整体符合低雷诺数翼上凸下凹的发展趋势。
图6 优化翼型与基础翼型对比
Fig.6 Comparison of optimized airfoil and the base airfoil
图7给出优化后翼型与基础翼型的表面压力系数对比。从图7中可以看出,前缘压力变化较大,横坐标从0.05到0.4左右前缘上表面压力得到较为明显的提升,0.2之前下表面压力有小幅度的回落。后缘变化为小幅扰动,变化并不明显。优化翼型与基础翼型相比,翼型上表面曲线后半段变化并不明显,与之相应的表面压力系数曲线两者近似拟合。翼型下表面曲线有小幅度的上凸,图7中从0.5到0.8左右相应长度压力系数得到小幅度平稳的提升。
图7 翼型表面压力系数分布对比
Fig.7 Comparison of pressure coefficient distribution between optimized airfoil and basic airfoil
基础翼型压力分布图与优化后翼型压力分布图如图8所示。可以看出优化后的翼型上表面曲线明显更弯,前缘受力点下移,翼型下表面后端出现压力集中。同时上表面所受压力和基础翼型相比明显向前缘移动,后缘所受压力并没有发生明显变化,综合来说优化后翼型上下压差增大,升阻比提高。
图8 优化翼型与基础翼型压力分布对比
Fig.8 Comparison of pressure distribution between optimized airfoil and basic airfoil
4.2 巡飞弹实例验证
为验证优化后翼型对于巡飞弹气动特性有相应的提升,现设计一款口径155 mm串列翼巡飞弹,弹体长度为930 mm,前后翼展长度分别为1 400 mm和1 200 mm,头部设计尺寸长度为247.5 mm,机身长600 mm。根据李永泽[16]提出的重叠式方案下装载口径与机翼弦长的关系式计算出机翼弦长为100 mm,巡飞弹原始模型如图9所示。为比较翼型优化对于巡飞弹巡飞性能的影响,其他部分不变,替换前后翼,弦长相同都为100 mm,翼型分别为E387翼型和优化翼型,经过气动仿真后进行比较分析,如图10所示。
图9 巡飞弹原始模型
Fig.9 Loitering munition original model
图10 分别采用基础翼型和优化翼型的巡飞弹设计
Fig.10 Loitering munition design with base airfoil and optimized airfoil
采用SCDM完成流场模型划分,如图11所示。其中圆柱形外边界为压力远场,巡飞弹表面为绝热无滑移壁面,采用Body of influence加密来精确捕捉流场变化。
图11 SCDM巡飞弹流场划分
Fig.11 Loitering munition fluid field decomposition in SCDM
巡飞弹设计飞行速度M=0.088,攻角α=4°,巡飞高度h=150 m,雷诺数Re=2×105。采用Fluentmeshing进行网格划分,通过Fluent进行气动特性分析导出升力系数及阻力系数,进行优化前后比较分析。
如图12所示,可以较为明显的看到在机身其他部分压力保持平稳的同时,采用优化翼型的巡飞弹机翼上表面压力更大,说明优化翼型巡飞弹上下翼压力差相较于采用基础翼型的巡飞弹更大,获得的升阻比更高,气动性能进一步提升,填补了巡飞弹低雷诺数下翼型的空白。
图12 优化翼型与基础翼型巡飞弹压力分布对比
Fig.12 Comparison of pressure distribution between loitering munition with optimized airfoil and loitering munition with basic airfoil
气动计算后,基础翼型的升力系数为0.606 00,阻力系数为0.071 92,升阻比为8.425;优化翼型的升力系数为0.663 76,阻力系数为0.074 09,升阻比为8.958,如表2所示。对比而言,采用优化翼型的巡飞弹气动性能够得到提升,升阻比提高了0.533,提升百分比为6.326%,升力系数提高了0.063 76,提升百分比为9.531%,证实了本文中设计的基于CST参数化方法的巡飞弹翼型优化设计对巡飞弹设计有实际参考意义。
表2 优化前后巡飞弹气动力系数值
Table 2 Loitering munition aerodynamic coefficients before and after optimization
巡飞弹模型ClCdCl/CdΔCl/Cd采用基础翼型0.606 000.071 928.425采用优化翼型0.663 760.074 098.9586.326
5 结论
1) 对低雷诺数E387翼型CST参数化方法进行研究,发现Bernstein 多项式阶数n>6时精度基本达到仿真要求,8阶以后,误差趋于稳定,可以达到较高的精度。
2) 对低雷诺数翼型E387优化设计后,获得的优化翼型的升阻特性在巡飞条件下得到了有效改善,升阻比提高了1.917 3,提高百分比为4.173%,升力系数提高了0.060 19,提高百分比为11.127%。
3) 采用优化翼型的巡飞弹气动性能得到进一步提升,升阻比提高了0.533,提高百分比为6.326%,升力系数提高了0.063 76,提高百分比为9.531%,满足巡飞弹对于特定翼型的需求,为低雷诺数巡飞弹提供理论依据和设计参考。
[1] MARK V.Performance analysis and design of loitering munitions:A comprehensive technical survey of recent developments[J].Defence Technology,2022,18(3):325-343.
[2] 朱志斌,刘强,白鹏.低雷诺数翼型层流分离现象大涡模拟方法[J].空气动力学学报,2019,37(6):915-923.
ZHU Zhibing,LIU Qiang,BAI Peng.Large eddy simulation of laminar separation in low Reynolds number airfoils[J].Acta Aerodynamica Sinica,2019,37(6):915-923.
[3] SELIG M,GUGLIELMO J,BROERN A,et al.Experiments on airfoils at low Reynolds numbers[C]//34th Aerospace Sciences Meeting and Exhibit.Reno:AAIA,2013.
[4] BREHMC,MACK S,GROSS A,et al.Investigations of an airfoil at low reynolds number conditions[C]//4th Flow Control Conference,2008.
[5] 廖炎平,刘莉,龙腾.几种翼型参数化方法研究[J].弹箭与制导学报,2011,31(3):160-164.
LIAO Yanping,LIU Li,LONG Teng.Study on several airfoil parameterization methods[J].Journal of Guidance,Projectiles,Rockets,Missiles andGuidance,2011,31(3):160-164.
[6] WALTER M,FRANCISCO J.P.A sequential quadratic programming algorithm using an incomplete solution of the subproblem[J].SIAM Journal on Optimization,2006,5(3):590-640.
[7] 钱瑞战,乔志德,宋文萍,等.基于N-S方程和序列二次规划的翼型优化设计[J].西北工业大学报,2003(5):523-527.
QIAN Ruizhan,QIAO Zhide,SHONG Wenping,et al.Optimal Design of airfoil based on N-S equation and sequential quadratic programming[J].Journal of Northwestern Polytechnical University,2003(5):523-527.
[8] KULFAN B,BUSSOLETTIJ."Fundamental" parameteric geometry representations for aircraft component shapes[C]//11th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference,2006.
[9] KULFAN B.A universal parametric geometry representation method-“CST”[C]//45th AIAA Aerospace Sciences Meeting and Exhibit,2007.
[10]关晓辉,李占科,宋笔锋.CST气动外形参数化方法研究[J].航空学报,2012,33(4):625-633.
GUAN Xiaohui,LI Zhanke,SHONG Bifeng.Research on parametrization method of CST aerodynamic profile[J].Chinese Journal of Aeronautics,2012,33(4):625-633.
[11]KUIFAN B.Universal parametric geometry representation method[J].Journal of aircraft,2008,45(1):142-158.
[12]CEZEM,HAYASHI M,VOLPE E.A study of the CST parameterization characteristics[C]//27th AIAA Applied Aerodynamics,2009.
[13]卜月鹏,宋文萍,韩忠华,等.基于CST参数化方法的翼型气动优化设计[J].西北工业大学学报,2013,3(5):829-836.
PU Yuepeng,SHONG Wengping,HAN Zhonghua,et al.Aerodynamic optimization design of airfoil based on CST parameterization method[J].Journal of Northwestern Polytechnical University,2013,31(5):829-836.
[14]RAFAEL P T.A numerical study on the unsteady aerodynamic characteristics of thick aerofoils with wavy leadingedges[D].Southampton:University of Southampton,2016.
[15]SELIG M S,MCGRANAHAN B D.Wind tunnel aerodynamic tests of six airfoils for use on small wind turbines[C]//42nd AIAA Aerospace Science Meeting and Exhibit.Reno:AIAA,2004:1-19.
[16]李永泽,孙传杰,卢永刚.巡飞器载重能力与装载火箭弹间匹配关系研究[J].飞行力学,2013,31(6):570-573.
LI Yongze,SUN Chuanjie J,LU Yonggang.Study on the matching relationship between the payload capacity and the loaded rocket of the aircraft cruiser[J].Flight Dynamics,2013,31(6):570-573.