近年来，随着能源紧缺和环境污染问题的日益突出，世界各国对零污染，大储量、高效率的氢能产生极大兴趣，各国也在大力推动氢能源飞机的研制[1-2]。与传统无人机相比，氢能无人机具有高效、无污染、长续航等优点[3]。然而只以氢能源为动力驱动的无人机，由于燃料电池本身的特性，需要采用非常规气动布局、大展弦比柔性机翼等技术来适配氢燃料电池额外带来的空间布局[4]，这样也会造成设计成本高、设计周期长和飞行性能下降等问题[5]。因此，针对绿色能源，长航时的设计要求，采用混合动力驱动的无人机因其成本低、能执行低空长航时任务成为无人机领域的研究热点。
翼型选型是长航时氢能源无人机机翼设计的首要环节，各国学者针对翼型优化问题进行了大量的研究。陈杰等[6]采用基于自组织特征映射(SOM)神经网络的超临界翼型设计方法，研究了超临界翼型设计问题。赵童等[7]考虑超临界机翼的当地后掠角和当地曲率影响，发展了面向三维机翼性能的翼型优化设计方法。刘远强等[8]等针对某型通用飞机设计状态，基于自由变形(FFD)技术参数化方法以及NSGA-Ⅱ多目标优化算法进行了高升力层流翼型的优化设计，所提出的方法可为通用航空飞机高升力层流翼型设计提供参考。赵欢等[9]通过对影响气动稳健优化设计效率的关键技术进行研究，发展了基于自适应前向-后向选择(AFBS)的稀疏多项式混沌重构方法，有效地解决了传统翼型设计方法难以满足高速高升力自然层流翼型设计兼顾高升力设计、自然层流设计以及超临界设计的难题。
气动特性是表征飞行器飞行性能的重要指标，常采用计算流体动力学(CFD)和其它优化算法相结合的方法对其进行优化分析。吴文华等[10]基于伴随算子的气动布局优化方法，开展全机状态下的机翼多参数、高精度优化设计取得了明显的优化效果。在国外，OKTAY et al[11]针对飞机机翼结构拓扑优化问题，结合CFD，提出了一套结构优化工具。MASOUD et al[12]通过风洞试验和CFD建模，研究运输机尾部进行降阻增升优化，结果表明CFD可作为阻力计算和飞机外形设计的可靠工具。BOELENS et al[13]的研究结果也表明CFD技术在飞机设计过程中起到不可替代的作用。YUE et al[14]利用CATIA软件设计了机载加油机的三面几何模型，并基于CFD技术对装载式加油机流场的气动特性进行了数值模拟，得到了压力、速度和喷流温度等流场的基本特性。SURESH et al[15]为了提高单翼布局的效率，提出了“C翼”布局，利用CFD技术比较了直翼和“C翼”的阻力、失速等特性。但目前针对氢能源长航时无人机的气动分析较少，其气动布局特性有待进一步研究。
本文以长航时氢能源混合动力无人机为研究对象，基于MH114翼型进行优化设计，得到满足研制需求的翼型，提出直机翼(HY0)和椭圆翼(HY1)两种气动布局，并对两种布局分别处于最大升阻比和失速两种状态下的流场进行数值仿真，分析流线、压力分布、气动参数等流场特征，以期评估该款无人机气动布局方案，为进一步优化设计和改进提供参考。
1 构型设计方案
长航时氢能源混合动力无人机主要性能指标如表1所示。为了满足设计要求，该款氢能源动力无人机选择整体强度大、稳定性好的双尾撑布局，前三点式起落架。固定翼选取大展弦比、上单翼的布局形式，基于MH114翼型对其进行改进优化得到综合性能更好的翼型；动力储氢装置(含燃料)质量约322.3 kg，容积为753.6 L，氢燃料重量约29.03 kg，总动力装置约416.72 kg，巡航螺旋桨功率约为10.7 kW，燃料电池电堆的输出功率约为14.7 kW，经过参数匹配，燃料型号巴拉德公司产 FCvelocity®-9SSL较符合；考虑具有不同平面形状的直机翼与椭圆翼两种布局方案，分别对其进行气动分析以获得性能更好的机翼形状。双尾撑尾端各布置一个垂尾，双垂尾和平尾呈“Π”型布局，减少机翼尾涡带来的不利干扰，为飞机提供足够的稳定性和操纵性。该无人机整体呈现出“廿”式布局，三维模型如图1所示，几何尺寸参数如表2所示。该无人机起飞速度v1=79.2 km/h，巡航速度v2=100 km/h，巡航高度1 km，考虑起飞、速度变化(蓄电池工作)等工况，续航阶段航时为30 h/90%；海平面的最大平飞速度为vmax=164.6 km/h，最小平飞速度为vmin=71.8 km/h，净升限为4.2 km.该无人机的任务剖面图见图2.其中，2点和3点分别为巡航开始点和结束点，v为飞行速度、H为飞行高度，CL为升力系数，设整个续航阶段V、H和CL均相等。
表1 无人机性能要求
Table 1 UAV performance requirements
主要参数性能要求最大起飞重量≤800 kg有效载荷≥150 kg最大航时≥30 h最大航程≥3 000 km巡航速度≥100 km/h巡航高度≥1 000 m
表2 无人机主要几何参数
Table 2 Geometric parameters of UAV
机型机长/m机高/m机翼面积/m2展弦比平均气动弦长/mHY0102.615.663200.861HY1102.617.416200.918
图1 两种无人机气动布局
Fig.1 Two pneumatic layout of UAV
图2 无人机任务剖面图
Fig.2 UAV mission profile
2 数值模拟
2.1 网格划分
考虑到无人机模型的结构对称性，为提高数值计算速率，建立半机模型并对其外部流场采用Pointwise进行网格划分。整个计算区域为半个立方体，区域的前部到头部机尾的距离为30L，上下面距离机身15L，展向距离机翼30L，其中L为机身长度，如图3(a)所示。网格的生成对无人机进行气动分析至关重要，网格质量直接影响求解效率和计算精度。对整机进行面网格划分时，机翼前后缘以及模型曲率较大区域做加密处理，在此基础之上生成体网格，其中整机区域附近网格细化且边界层体网格设置为100层，并保证第一层网格无量纲距离y+≈1，以满足湍流模型的计算要求。无人机网格分布如图3所示，经过网格无关性验证优化后网格总数约为690万。
图3 长航时无人机网格划分示意图
Fig.3 Schematic diagram of grid division of long-endurance UAV
2.2 控制方程
控制方程采用稳态的时均N-S方程，包括连续方程、动量方程和能量方程，具体可表示为如下通用形式：
(1)
式中：ρ为密度，ν为速度矢量，φ为待求变量，包括速度、温度等Γφ为广义扩散系数，Sφ为广义源项。广义扩散系数和广义源项的具体表达式由待求变量决定。
2.3 数值方法与边界条件
本文采用ANSYS CFX软件对流场进行数值计算，计算中使用整场离散耦合求解的策略，应用有限体积法离散控制方程组，对流项和扩散项的离散采用2阶迎风格式。采用双参数SSTk-ω(Shear Stress Transfer)湍流模型求解时均N-S方程。SSTk-ω模型综合了k-ω模型在近壁区计算的优点和k-ε模型在远场计算的优点，同时增加了横向耗散导数项，在湍流黏度定义中考虑了湍流剪切应力的输运过程，适用更广，更加适用于带逆压梯度的流动、翼型气动分析等问题。SSTk-ω湍流方程为
(2)
(3)
式中：Gk是由层流速度而产生的湍流动能；Gω是由ω方程产生，Γk和Γω表明k和ω的扩散率；Yk和Yω是由于扩散产生的湍流；Sk和Sω为表面张力。
来流风速为100 km/h，高度为1 000 m，外部大气压为89 875.8 Pa，静温为281.65 K，侧滑角为0°，以经验公式选取最大升阻比迎角α=4°和失速迎角α=12°两种工况进行数值模拟。流场中气体为理想气体，入口、出口和其他远场均选择自由流边界条件，对称面设置对称边界条件，飞机表面选择无滑移壁面边界条件。
3 结果与讨论
本文首先对翼型进行优化，然后建立模型对其进行流场数值模拟，得到两种气动布局下的速度场、流线、压力和马赫数分布等特征，重点对翼型优化、机翼表面流场和压力分布以及飞机气动参数进行分析。
3.1 翼型优化
基于MH114基础翼型，联合采用MATLAB和XFOIL软件，使用Hicks-Henne型函数法对翼型进行参数化，以雷诺数为1 820 000、迎角4°时的升阻比为优化目标，使用遗传算法对基础翼型进行优化，适应度函数利用XFOIL作函数器求解升阻比，设计变量为Hicks-Henne型函数的12个系数。得到了优化前后的翼型曲线对比、升力系数对比、阻力系数对比、升阻比对比和俯仰力矩系数对比，分别如图4(a)-(e)所示。从图中可以看出迎角为4°时，升力系数由1.306 7提高到了1.468 5，提高了12.38%；阻力系数由0.007 15增加到0.007 19，增加了0.42%；俯仰力矩系数绝对值由0.191 8增加到0.227，增加了18.35%；升阻比由180.336提高到205.644，提高了14.03%；因此，迎角为4°时优化后的翼型性能有了很大的改善。在整个迎角变化过程中，升阻比在-8°到7°之间都高于基础翼型，阻力系数在-8°到8°之间与基础翼型基本相同，所以优化后在无人机平飞状态时的迎角范围内翼型阻力基本没有变化，而升阻比有了显著的提高。
3.2 机翼表面流场分析
机翼是飞机产生升力的主要部件，无人机整体性能的优劣往往与机翼表面气流的附着情况有关。为了直观地分析两种构型机翼在不同迎角状态下的流场流动情况，图4和图5分别给出了在最大升阻比状态(迎角α=4°)和失速状态(迎角α=12°)下两种构型机翼表面的局部外流场气流分布，为了更好地探究不同迎角下流动分离的机理，在距离无人机中性面3 m位置，分别取机翼纵截面外流场流线进行分析。从图5中可以看出α=4°时，两种构型无人机外流场表面流线依然为附着流动；从图6中可以看出随着飞机抬头，迎角不断增大，α=12°时下后缘流线均出现了大面积的分离，并有不断向前缘发展的趋势，飞机开始出现失速。
图4 优化前(MH114)、后(optim)翼型气动参数对比
Fig.4 Comparison of aerodynamic parameters of airfoils before and after optimization
图5α=4°时，两种构型无人机表面流场对比
Fig.5 Airfoil surface flow field comparison whenα=4°
图6α=12°时，两种构型无人机表面流场对比
Fig.6 Airfoil surface flow field comparison whenα=12°
3.3 机翼表面的压力分布
图7展示了α=4°和α=12°下，两种构型无人机表面压力分布云图。从图中可以看出，两构型的整体压力布局较合理，最大压力普遍分布在机翼的前缘下表面，最小压力也集中在机翼上表面的前缘位置，这就使全机压差的最大值主要集中在靠近机翼前缘的大部分区域，尽可能的为无人机贡献更多的升力。同时，α=4°，机翼上表面展向压力分布分界线近似直线分布，过渡十分平缓，说明该类压力分布主要由翼型形状决定，机翼表面气流未发生流动分离，并且椭圆翼构型无人机的升力趋于翼根处集中，分布更合理。当飞机处于失速状态α=12°时，机翼上表面压力分布界线分布不再均匀，说明在大迎角飞行时，机翼表面的展向流动加剧，机翼表面已经开始了流动分离；进一步对比两种机型的压力分布，直机翼构型的波浪状压力界限位置相比于椭圆翼构型更加靠前，椭圆翼构型的失速特性和升力分布要优于直机翼构型。
3.4 气动参数
图7给出直机翼和椭圆翼构型无人机的基本气动力系数CFD分析结果，结果显示，随着迎角的增加，两构型的升力系数先近似呈线性增长，直形翼构型无人机迎角α=10°左右，升力系数到达最大值，开始发生失速；椭圆翼构型无人机迎角α=12°左右，进入失速状态，也再次证明椭圆翼构型的失速特性要稍优于直机翼构型；从力矩特性来看，两构型无人机都满足静稳定状态，但椭圆翼构型的力矩系数较小，力矩特性要优于梯形翼构型。由于两构型的整体尺寸、结构布局相差不大，所以升力和阻力数值上相差不大，但升阻比都到达20.3左右，整体气动性能满足本次长航时氢能源混合动力无人机的设计要求。综合来看，椭圆翼构型从稳定性、失速特性、升力分布以及强度上要优于直机翼。
图7 机背与机腹压力分布比较
Fig.7 Comparison of pressure distribution between back and belly
图8 气动参数对比
Fig.8 Comparison of pneumatic parameters
通过气动分析结果可知，HY1(椭圆翼)飞机巡航推重比为0.049 2，功重比为1.37 m/s，巡航时螺旋桨功率为10.7 kW，飞机以100 km/h的速度巡航30 h，考虑工况变化，实用航时为30/90%.计算巡航阶段所需能量为1 284 000 kJ.氢气的热值为1.4×105kJ/kg，氢燃料电池发电效率约为50%，估算螺旋桨效率为85%，电机电调效率为90%，直流变换器DC/DC效率为95%.则计算所需氢燃料质量为25.24 kg.另外还需要考虑起飞、速度变化(蓄电池工作)等工况，要适当考虑燃料裕度，一般增加10%～20%余量(取15%).因此所需氢气质量为29.03 kg.以70 MPa、质量储氢密度9%来计算氢容器的质量，则氢气罐(含燃料)质量约为322.6 kg，容积为753.6 L，总的动力装置重量为416.72 kg.加上飞机结构、飞控及通讯系统的重量，无人机总重量为725.28 kg.无人机设计满足航时不低于30 h的情况下，也同时满足最大起飞重量，巡航高度以及航程设计要求。HY0(直翼)验证如上。因此，无人机整体气动性能满足本次长航时氢能源混合动力无人机的设计要求。
4 结论
本文针对两种满足性能指标的氢能源混合动力固定翼无人机气动布局进行气动分析与优化，对两种气动布局方案进行评估，为进一步优化提供了方向，通过本文研究可以得到如下结论：
1) 通过联合MATLAB和XFOIL对MH114翼型进行优化，优化后的升阻比由180.336提高到了205.644提高了14.03%，翼型的气动性能相较基础翼型有了很大的提高，更适合本项目的需求。
2) 最大升阻比状态(迎角α=4°)，两构型飞机翼身连接处后缘均发生不同程度流动分离；失速状态(迎角α=12°)下，机翼上表面压力分布不均匀，机翼后缘出现分离涡，但HY1方案气流依附更加紧密。
3) 两种构型无人机的压力分布较合理，最大压差主要集中在机翼前缘；最大升阻比状态(迎角α=4°)，HY1升力趋于翼根处集中，分布更好；失速状态(迎角α=12°)下，HY1构型的波浪状压力界线位置更靠后，说明HY1构型的失速要优于HY0.
4) 从力矩特性来看，两构型无人机都满足静稳定状态，但HY1构型的力矩系数较小，稳定性要优于HY0构型。
[1] 赖耀胜,李龙.氢能飞机发展现状分析[J].航空动力,2021(6):37-40.
LAI Y S,LI L.Hydrogen powered aircraft[J].Aerospace Power,2021(6):37-40.
[2] 徐东,刘岩,李志勇,等.氢能开发利用经济性研究综述[J].油气与新能源,2021,33(2):50-56.
XU D,LIU Y,LI Z Y,et al.A cost-efficiency review of hydrogen energy exploitation[J].Petroleum And New Energy,2021,33(2):50-56.
[3] 田海平,苏健,肖泽琦,等.氢动力固定翼无人机气动布局的气动特性评估[J].太原理工大学学报,2021,52(4):607-613.
TIAN H P,SU J,XIAO Z Q,et al.Aerodynamic layout evaluation of hydrogen-powered fixed-wing unmanned aerial vehicle(UAV) by numerical simulation[J].Journal of Taiyuan University of Technology,2021,52(4):607-613.
[4] 李延平.太阳能/氢能混合动力小型无人机总体设计[D].北京：北京理工大学,2014.
[5] 李延平,刘莉.太阳能/氢能混合动力无人机及关键技术[J].飞航导弹,2014(7):39-45.
[6] 陈杰,孙刚.基于SOM神经网络的超临界翼型设计[J].力学季刊,2011,32(3):411-417.
CHEN J,SUN G.Supercritical airfoil design based on SOM neural network[J].Chinese Quarterly of Mechanics,2011,32(3):411-417.
[7] 赵童,张宇飞,陈海昕,等.面向三维机翼性能的超临界翼型优化设计方法[J].中国科学:物理学力学天文学,2015,45(10):89-101.
ZHAO T,ZHANG Y F,CHEN H X,et al.Aerodynamic optimization method of supercritical airfoil geared to the performance of swept and tapered wing[J].Sci Sin-Phys Mech Astron,2015,45(10):89-101.
[8] 刘远强,李天,白俊强,等.通用飞机高升力层流翼型优化设计研究[J].西北工业大学学报,2017,35(2):339-347.
LIU Y Q,LI T,BAI J Q,et al.Optimization design of high-lift laminar airfoil for general aircraft[J].Journal of Northwestern Polytechnical University,2017,35(2):339-347.
[9] 赵欢,高正红,夏露.高速自然层流翼型高效气动稳健优化设计方法[J].航空学报,2022,43(1):271-288.
ZHAO H,GAO Z H,XIA L.Efficient robust aerodynamic design optimization method for high-speed NLF airfoil[J].Acta Aeronautica et Astronautica Sinica,2022,43(1):271-288.
[10] 吴文华,范召林,陈德华,等.基于伴随算子的大飞机气动布局精细优化设计[J].空气动力学学报,2012,30(6):719-724.
WU W H,FAN Z L,CHEN D H,et al.Adjoint based on high precise aerodynamic shape optimization for transonic civil aircraft[J].Acta Aerodynamica Sinica,2012,30(6):719-724.
[11] OKTAY E,AKAY H U,MERTTOPCUOGLU O.Parallelized structural topology optimization and CFD coupling for design of aircraft wing structures[J].Merttopcuoglu Computers & FluidsVolume,2011,49(1):141-145.
[12] MASOUD MIRZAEI,MOHAMMAD HADI KARIMI,MOHAMMAD ALI VAZIRI.An investigation of a tactical cargo aircraft aft body drag reduction based on CFD analysis and wind tunnel tests[J].Aerospace Science and Technology,2012,23(1):263-269.
[13] BOELENS O J.CFD analysis of the flow around the X-31 aircraft at high angle of attack[J].Aerospace Science and Technology,2012,20(1):38-51.
[14] YUE K Z,CHENG L L,ZHANG T Y,et al.Numerical simulation of the aerodynamic influence of an aircraft on the hose-refueling system during aerial refueling operations[J].Aerospace Science and Technology,2015,49:34-40.
[15] SURESH C,RAMESH K,PARAMAGURU V.Aerodynamic performance analysis of a non-planar C-wing using CFD[J].Aerospace Science and Technology,2015,40:56-61.