临近空间飞行器一般指其主要活动区域在20~100 km，此范围近年来成为各国开展高速飞行器研究的主要领域。飞行器气动外形设计的目的是设计最合理的气动外形，使得飞行器在给定的约束条件下获得最优良的气动性能^[1-2]。
在临近空间范围内，为兼顾装填比及升阻特性需求，考虑工程化和实用性，飞行器以升力体外形为主，在布局设计中，需考虑诸多限制和约束因素，对飞行器的选型和优化常耗费大量时间，依靠个人设计经验调节参数优选方案的传统气动布局设计方式面临巨大的挑战，为提高设计效率，基于参数化外形的气动布局优化设计技术逐渐发展起来，为复杂气动外形设计问题提供了有效的解决途径^[3-7]。目前，国内研究方面，西北工业大学和南京航空航天大学主要以低速航空飞行器为主^[8-11]，北京航空航天大学、中国空气动力研究与发展中心以及中国航天空气动力技术研究院聚焦于高速航天飞行器^[12-15]，针对气动布局的优化设计主要集中于机翼或一体化体身等部件的外形参数优化，参数较少，并且未充分考虑到工程化的约束条件。为推进飞行器的设计技术研究，本文探索了一套优化设计方法，以主要总体性能指标为约束条件和优化目标，通过几何参数的相关性分析及排序，针对关键参数的优化迭代，形成符合总体需求的气动外形，为布局设计的快速优化迭代提供新的技术手段。
1 总体性能相关性分析 在飞行器气动布局设计中，为追求升阻特性及飞行性能的提升，广泛应用了各种优化设计技术。这个过程离不开合理的飞行器参数化模型和适当精度的气动特性评估手段。前者需要一种恰当的参数化几何建模方法，用尽可能少且物理意义明确的设计参数来描述飞行器的几何外形，当设计参数变化时，需要方便快捷地得到新的飞行器外形；后者则需要一种兼顾效率和精度的气动特性计算方法，可以高效、快速地估算飞行器的气动特性，并且在宽速域范围内满足一定精度要求。
1.1 几何外形参数化建模 在进行气动外形优化计算之前，需要将相关的几何外形使用有限数量的参数表示出来，外形参数化是气动外形设计和优化的第一步，并且对优化过程有着十分重要的影响。由于搜索设计空间的计算量随着设计参数的增加而呈指数增长，参数化过程应使用尽可能少的设计参数来表示几何外形，以降低优化过程的计算量。
高速飞行器目前以翼身组合的面对称升力体为主，考虑装填比和升阻比需求，构造了一个具有普适性的初始几何外形(见图 1)。其主体为双锥升力体布局，主要由球形端头、双锥拱形身部、主翼、升降舵和方向舵等部分组成。图中：BS、BA和BB分别分别表示相对应的截面编号；R_h为飞行器的钝化半径；L为飞行器总长；L₁表示第一锥段的长度；L_W为机翼的展长。第一机翼段包络线为直线，第二机翼段包络线为卡门曲线。机翼特征横截面有3个，分别为WS、WA和WB横截面。为了便于优化设计，将两段机翼合并为一段，并将在xOz平面内的展长包络线设计为二次曲线。由于该飞行器表面外形较复杂，为了便于利用面元法进行气动特性计算，参数化建模的基本原则是：在不影响整体形状的前提下，忽略不便于建模的次要细节，尽量完全复现飞行器的真实外形，同时提高参数化建模的快速性与高效性。


根据对初始外形特征的简化调整和描述提取，归纳出整个飞行器在几何参数化建模时所需的49个设计参数，根据几何设计参数在量纲上的差别和影响程度的大小，分为3类，即长度量纲类(单位：mm)、角度量纲类(单位：(°))和系数无量纲类，部分几何参数见表 1。
表 1 参数化建模的主要几何参数 Table 1 Main geometrical parameters of parametric modelling


表选项






在总体方案论证过程中，气动设计不能单纯追求高速下的高升阻比指标，在飞行器参数化几何外形的基础上，为有效分析参数影响度，基于总体指标的合理性以及考虑装填比、防隔热和结构等各专业因素(见表 1)，对参数进行了取值约束。
1.2 快速工程算法评估 在对飞行器的基本气动特性预示过程中，采用面元法划分复杂外形的表面，通过将平板的摩阻加入形状因子和可压缩因子来获得飞行器的摩阻。机身迎风面采用Dahlem-Buck公式，机身背风面采用Prandtl-Meyer公式，机翼迎风面采用切楔/切锥法，机翼背风面采用膨胀波方法。
采用快速工程算法对参数化外形的升阻特性及舵效进行了评估，得到一套典型升力体的气动特性数据库，与风洞试验的结果进行了比对，图 2摘取了Ma=10.5时快速工程算法得到的计算数据和风洞试验数据。图中：α为攻角；β为侧滑角。可见，方向舵δ_n偏转对轴向力系数C_A及法向力系数C_N的影响不大，由比对数据分析，轴向力系数的偏差为2%~5%，法向力系数的偏差最大为9%，可见轴、法向力的吻合度满足评估需求。考虑到在飞行器气动数据库建设过程中，除对效率的高要求以外，核心考核指标为：在偏差范围内，计算数据在对气动力及力矩系数变化趋势的跟踪上能保持一致。由此评判，所采用的工程算法可开展批量数据的快速生成，搭建起飞行器的气动特性数据库。


1.3 几何外形参数相关性排序 几何设计参数变化将直接改变飞行器的几何外形，从而间接改变飞行器的气动特性，选择最关注的几个总体性能指标进行分析，分别为阻力系数f₁=C_D、升阻比f₂=κ、容积率f₃=η、纵向静稳定裕度f₄=X_cp-X_g以及操纵特性f₅=m_z^α/m_z^δ_z，X_cp为压心位置，X_g为重心位置，m_z^α和m_z^δ_z分别为俯仰力矩系数对攻角的偏导数和俯仰力矩系数对舵偏角的偏导数。f_GEO表示与几何外形相关的性能指标参数，f_FLY表示与飞行性能相关的性能指标参数，H表示飞行高度。因此，通过计算可以获得总体性能指标与设计参数之间的隐式关系，如下：

x∈Rⁿ⁼⁴⁹为几何设计参数变量，包含所有几何设计参数。
参数影响度的大小用灵敏度函数来表征。基于一阶灵敏度函数的参数影响度分析方法(见表 2)，可有效表征各个设计参数对系统响应的影响，进行敏感度分析时，马赫数取值范围为5 < Ma < 10.5。
表 2 一阶灵敏度函数 Table 2 First-order sensitivity function


表选项






表 2中：f为系统特性变量；x为参数变量。当参数变量x各分量的量纲一致时，灵敏度函数计算采用一阶标准形式；当参数变量x的增量Δx以百分比形式给定时，灵敏度函数计算采用一阶百分比形式；当参数变量x各分量的量纲不一致时，灵敏度函数计算采用一阶对数形式。将飞行器特性变量y_i=f_i(x)(i=1, 2, …, n)对参数变量x=[x₁, x₂, …, x_r]在基准状态x=x₀处的一阶偏导数作为评判变量x中各元素对飞行器特性影响程度大小的依据。在基准状态x=x₀处给定增量Δx，运用有限差分原理，求得一阶灵敏度函数S₁为

给定基准状态x=x₀和参数增量Δx，采用一阶对数形式求取灵敏度函数，对总体性能与气动参数进行相关度分析，对影响总体性能指标的关键几何参数进行了排序，其中，影响飞行器升阻比的几何参数如表 3所示。
表 3 参数敏感性排序 Table 3 Rank of sensitive parameters on each factor


表选项






为使分析结果更具有可比性和合理性，按照长度量纲、角度量纲和系数无量纲系统进行分类评估，由结果发现，高升阻比布局的设计方向应聚焦一锥截面及主翼几何外形，通过控制端头半径、一锥锥角及空气舵几何外形进行局部减阻，优化主翼形线及安装角达到增升目的；针对飞行器稳定性的几何参数主要关注一锥锥角、主翼几何外形及长度系数，体身设计完成后，通过调整空气舵展长及舵位置系数来调节静稳定裕度。
2 面对称飞行器优化平台设计 飞行器布局优化设计是一个多学科优化问题，CFD技术与数值优化方法结合在一起，实现基于CFD的飞行器气动外形优化设计^[16-19]，建立计算精度高、适应性强、使用方便、人工干预少和稳定可靠的设计优化环境或框架是未来发展趋势。
优化主要思路为：基于飞行器相关性分析的几何参数排序结果，运用Kriging代理模型以及DIRECT方法进行气动优化设计，过程中引入动态网格技术，包括表面网格运动技术和体网格运动技术，实现在优化器得到气动外形后CFD计算网格的自动处理，最终实现优化设计自动化。
2.1 Kriging代理模型 当目标函数评估比较耗时时，在优化过程中用一个合理的数学模型来代替真实的目标函数评估，这个数学模型称为代理模型^[20-23]。Kriging方法是1951年由Krige^[24]提出的一种估计方差最小的无偏估计模型，其通过相关函数的作用，具有局部估计的特点，对非线性模型有很好的逼近能力，充分考虑了样本点之间的距离和样本点的团聚，在估值的过程当中能够很好地体现变量的空间分布情况与连续性。
在构造Kriging代理模型之前，采用拉丁超立方体抽样对样本空间进行采样。针对样本均匀性评估准则问题，采用Hickernell^[25]提出的L2中心差异方法对样本分布均匀性进行分析。过程中针对多种复杂函数进行了研究，验证了Kriging代理模型和样本生成方法模拟复杂非线性函数的能力和精度。
2.2 数值优化方法 气动外形优化设计过程的本质为：在约束条件下，在设计参数空间中寻找气动参数最值，在数学上即为约束条件下的数值优化问题。研究过程中采用DIRECT方法^[26]，其是一种直接全局优化方法，优化空间被剖分为超矩形(hyper rectangle)，在超矩形中心计算目标函数值，迭代过程中，优化算法选取适当的超矩形作为样本对象进行进一步剖分。
基于DIRECT方法的多维空间寻优方法具有编写简单、计算效率高、全局寻优能力强及鲁棒性强的特点，可通过结合基于Kriging方法的代理模型方法实现复杂飞行器气动布局要素优化设计目的。
2.3 优化平台集成及应用测试 由基于CFD的数值气动优化设计过程可以看出，在优化设计前需要完成气动外形的参数化工作，并且在优化中或者样本点建立过程中随着气动外形参数变化流场计算网格需要调整，而二者的鲁棒性直接影响优化设计方法的效率和鲁棒性。动网格采用基于Delaunay图映射和弹簧准则的网格变形方法，参数化和动网格方法是影响方法是否具有工程应用价值的重要因素，鲁棒性强和自动化程度高的参数方法和动网格方法可以有效减少在气动优化设计中的人工干预工作量，使得优化方法真正能够在工程问题中得到应用。
基于CFD和数值优化方法的气动外形优化设计技术(见图 3)在实际工程设计应用中还存计算效率、寻优全局性、优化过程自动化程度及鲁棒性等问题。其中，计算效率和寻优全局性问题主要与CFD方法和数值优化方法有关；优化过程自动化程度和鲁棒性问题主要与气动外形参数化优化和自动动态网格处理方法有关。针对这些问题采用了如下技术途径：


1)采用基于Kriging方法的代理模型来描述设计空间上的气动特性分布，根据适当的样本点生成气动特性代理模型后，在优化设计过程中，新外形的气动特性通过气动代理模型进行评估，计算效率会大大提高。
2)在设计空间上寻优采用基于DIRECT方法的数值优化方法，该寻优方法不但具有全局性，而且具有较高的计算效率和鲁棒性；在进行多目标优化设计时，采用加权方法进行处理；约束以罚函数的形式出现在优化方法中。
3)气动外形参数在样本生成时采用基于Catia二次开发的方法。
4)动态网格处理采用基于Gridgen脚本开发的方法。
为测试优化设计方法及设计平台，采用该方法针对初始外形开展气动外形优化设计研究。将容积率、静稳定度以及空气舵舵面铰链力矩作为约束条件，以升阻比为优化目标，马赫数取值范围为5 < Ma < 10.5，在保证空气舵舵效的前提下，尽可能提高飞行器升阻比，并在滑翔段范围内保持升阻比的稳定性。经设计参数与设计目标的相关性分析及排序，按照影响度大小从49个几何参数中选取26个参数作为优化设计参数，其他外形参数在优化设计中保持不变。优化设计参数给定了设计变量的取值范围，即设计空间，生成260个样本点，经判断发现有16个无法生成外形，因此有效样本为244个。调用Catia的二次开发程序，生成每个样本点的外形，并输出相应的IGES文件，通过调用Gridgen脚本文件生成每个样本的网格，图 4(a)为其中一个样本网格示意图。经检查，244个网格中有8个网格存在负体积，通过人工调整可以使这些网格中的负体积单元消失。由于计算网格为块结构网格，其拓扑相同，因此在气动特性计算时，首先对其中一个样本点进行计算，然后以该样本流场作为初始流场进行其他样本气动特性的计算，由于初始流场良好，可以用很少求解步数获得收敛的结果，计算量大大减少。图 4(b)为模拟得到其中一个样本点的流场结果。


通过多轮优化得到了气动外形，图 5所示为优化得到的外形与初始外形对比。可以看出，在考虑到热防护以及装填空间约束的前提下，优化外形的一锥长度增大，半锥角减小，机翼形线后掠角增大，机身底面更加扁平，外形趋向于改善波阻的特点。


图 6给出了优化前后外形在设计状态(H=37 km，Ma=8，α=8°)沿对称面和翼舵干扰处的流场结构。可见，沿流向飞行器头部出现一道典型的弓形脱体激波，使飞行器头部出现局部高压，构成飞行器主要的阻力来源，通过针对半锥角及一锥长度的优化设计，头部波阻明显减小；主翼前缘也出现一道压缩波，翼后缘的膨胀波和舵前缘的压缩波有一定程度的交汇，说明翼舵流场出现一定程度干扰，通过针对主翼翼型、翼舵外形及空间布放位置的优化，在翼面波阻减小的同时，主翼背风面分离区也有所缩小。


3 风洞试验验证 为评估优化效果，开展气动特性验证，选取?300 mm低密度风洞进行优化前后2种外形的试验验证，考核优化后气动性能。
3.1 试验吹风状态 以来流黏性干扰系数为相似参数，试验模型采用7%缩比模型，参考面积0.000 9 m²，参考长度0.25 m。采用内式六分量测力天平，试验典型状态, 包括模拟高度H、总温T₀、总压P₀以及单位雷诺数Re_∞L如表 4所示。
表 4 典型马赫数状态 Table 4 Parameters of typical Mach number


表选项






3.2 试验结果分析 图 7为优化前后外形试验状态下的纹影图。对比图 6的数值流场，头部主激波的模拟吻合度很好，优化前后外形比较而言，优化后外形激波强度减弱。同时随着模型攻角的增大，主激波上下呈现不对称性，渐强的迎风面激波向模型腹部贴近。攻角增加到8°时，模型头部激波与迎风面方向舵产生激波干扰，方向舵前缘压缩波将主激波外推，当攻角达到16°时，模型头部主激波已打在方向舵前缘，这对舵面的热防护环境十分不利，在设计飞行轨迹时，需要考虑尽量避免大攻角飞行。


表 5给出了优化前后外形的气动性能对比。可以看出，优化得到的外形在设计点轴向力系数减小了32.6%，法向力系数减小了10.3%，针对优化目标的设计，升阻比提高了16.9%。可见，所采取的优化设计方法满足了总体设计要求，设计效果明显。
表 5 优化前后气动特性(Ma=8，α=8°) Table 5 Aerodynamic characteristics before and after optimization (Ma=8, α=8°)


表选项






图 8给出了典型滑翔段马赫数下优化前后外形气动系数的对比(实线为初始外形，虚线为优化外形)。可见，优化外形阻力系数均明显低于初始外形，阻力系数的减小量随马赫数增加而增大，量值在30%~50%之间；升力系数略低于初始外形，量值在20%~40%之间。结合图 5外形变化分析，端头及一锥的尺寸调整使得波阻迅速减小，同时主翼形线的改变在减小阻力的同时，由于升力面积的缩减，升力也同步减小；升阻比在滑翔攻角(α=6°~10°)增量明显，优化外形升阻比最大提升19%~26%，最大升阻比为3.78(Ma=5，α=6°)。由于升力面的调整，最大升阻比攻角由初始外形的8°前移到优化外形的6°，有利于飞行器的小攻角低阻飞行；初始外形纵向压心位置在0.54~0.59之间，随攻角增大先前移再后移，在滑翔攻角随马赫数增加变化范围4%，优化后纵向压心位置后移5%，使得飞行器纵向静稳定性增强，且随攻角增大线性后移，在滑翔攻角随马赫数增加变化范围2%，这一点尤为重要，对于大空域，宽速域飞行器而言，压心变化范围小会大大减轻控制系统对舵面偏转角度的需求，并减小配平状态的升阻比损失。


4 结论 为建立一套相对完整的高速升力体气动布局选型及优化方法，开展了以下工作：
1)采用快速工程算法构建参数化外形的气动数据库，利用相关性分析方法，通过无量纲化、分类别处理，获得几何参数的敏感度排序。
2)综合考虑敏感度和计算量，以优化目标和总体约束条件为基础，提取关键几何参数，通过代理模型和优化方法，获得符合要求的优化外形。
3)通过数值评估校核确定优化外形的气动性能，并通过风洞试验比对优化前后外形的气动特性，验证优化平台的有效性。
在国内相关研究的基础上，升力体外形布局优化方面有了一定技术进步，主要体现在：
1)考虑了一类普适面对称升力体的全机布局优化，几何参数众多，外形相对复杂。
2)针对此类飞行器，基于一阶灵敏度函数的参数影响度分析方法，提取出了关键几何参数的敏感度排序。
3)在容积率以及操稳特性等总体约束条件的前提下，建立起了基于Kriging代理模型以及DIRECT方法的布局优化设计流程和平台设计思路，为优化技术在高速飞行器选型设计上的实用性进行了有益探索。

参考文献