基于连续伴随方法的高超声速飞行器高精度气动优化

高超声速飞行器具有实现全球快速到达和低成本进入空间的潜力，目前各主要航天大国都正在开展对高超声速飞行器的深入研究。高超声速飞行环境下，飞行器受到气动力、热载荷严酷，高超声速飞行器机体和推进系统的设计面临推阻平衡、升重平衡等难题，机体/推进一体化设计方法是解决这些难题的有效方法[1-2]。目前的一体化设计方法主要是基于无黏理论，高超声速来流条件下黏性对于飞行器气动性能有重要影响，实际应用中还要权衡几何修型导致的性能恶化[3]。气动设计作为概念设计的起始点，性能优良的气动外形有利于整个总体设计方案的封闭性和可实现性。在考虑黏性的情况下，针对高超声速飞行器开展高精度气动优化可以在一体化设计的基础上进一步提高飞行器的气动性能。

高精度气动优化需要借助精细化CFD(Computational Fluid Dynamics)工具开展气动性能预测，优化过程通常采用进化策略或梯度策略[4]。进化策略基于遗传算法、粒子群算法等优化算法对整个设计空间进行全局寻优，对于多设计变量优化问题，庞大的CFD计算开销是面临的主要困难，代理模型在一定程度上可以降低CFD计算开销，但仍存在精度和泛化性问题[5-6]。梯度策略根据目标函数梯度信息对设计空间进行局部寻优，伴随方法具有计算量与设计变量数目无关的优势，是气动优化中高效获取目标函数梯度的方法，已成熟应用于亚、跨和低超声速气动优化问题[7-9]。但对于高超声速流动，伴随方程推导需要考虑黏性输运系数变分的贡献；同时高超声速流动中存在的强激波间断对伴随方程对流项离散格式的精度和鲁棒性提出了更高的要求。

针对伴随方法在高超声速优化问题中的应用，国内外研究人员进行了一系列探索。Copeland等在考虑高温化学非平衡效应的情况下推导了连续伴随方程，基于Steger-Warming格式对伴随方程对流项进行离散，在高马赫数来流下针对二维翼型对连续伴随方法的计算精准度进行了验证[10]；Kline等采用连续伴随方法针对高超声速进气道开展了优化研究，但忽略黏性输运系数变分对伴随方程的贡献，伴随方程对流项离散采用JST(Jameson-Schmidt-Turkel)中心差分格式[11]。宋红超等基于离散伴随方法对高超声速单边膨胀喷管开展了优化研究[12]；高昌等基于连续伴随方法对二维高超声速进气道开展了优化研究[13]。目前对于伴随方法在高超声速优化问题中的应用还处于初步探索阶段，已有的研究工作主要针对二维或简单流场，对伴随方法在三维复杂高超声速流动中的应用缺乏深入研究和验证。

为进一步探索伴随方法在高超声速优化问题中的应用，建立高超声速飞行器高精度气动优化方法，本文基于Navier-Stokes(NS)方程推导了连续伴随方程和气动力目标函数对应的边界条件和灵敏度公式，考虑层流黏性输运系数对伴随方程的贡献，探究基于二阶熵修正Roe格式的伴随方程对流项离散格式，构建适用于高超声速流动的连续伴随求解器；基于连续伴随求解器构建针对高超声速飞行器的高精度气动优化框架；通过二维翼型和Sanger飞行器机翼优化算例验证连续伴随求解器的计算精度，以Sanger飞行器升阻比为优化目标开展基于连续伴随方法的高超声速高精度气动优化方法的应用研究。

1 优化方法

1.1 连续伴随方程

伴随方法的基本思想是将流动控制方程作为优化问题的约束条件，对于气动力优化问题，目标函数可以表示为统一形式：

(1)

式中：J为伴随方程目标函数；Sw为壁面边界；j为伴随方程目标函数积分项；p为压力；σ为应力；d为气动力目标函数方向向量；I为单位矩阵；n为壁面法向量；dS为积分面元。

对于不同的气动力目标函数：

(2)

式中：d为气动力目标函数方向向量；α为攻角；β为侧滑角；C∞=ρ∞Aref/2为气动力无量纲因子，其中ρ∞为来流密度，u∞为来流速度，Aref为参考面积；CD和CL分别为阻力系数和升力系数。将NS方程作为约束条件：

(3)

式中：R为流动控制方程残差；U为流动控制方程守恒变量；Fc为对流项；Fv为黏性项；μ为黏性输运系数。通过Lagrange乘子Ψ将NS方程引入目标函数：

(4)

式中：D为计算域；dV为积分体积微元。对式(4)进行变分运算：

(5)

式中：δ为变分运算符；δSw为壁面变分。目标函数的变分项根据曲面变分公式可以简化为

(6)

对于式(5)中NS方程的变分项，由于黏性项中包含流场变量的梯度，为了简化推导过程，将NS方程视为流场变量和流场变量梯度的函数：

(7)

对式(7)进行变分运算，根据求导交换法则：

(8)

式中：各Jacobian矩阵可表示为

(9)

式中：Ac为无黏通量Jacobian矩阵；Av为黏性通量Jacobian矩阵；Aμ为黏性输运系数变分项对伴随方程的贡献；D为黏性项关于守恒变量梯度的Jacobian矩阵。对于高超声速流动，黏性输运系数对近壁面流动起主导作用，对壁面压力和摩阻分布有显著影响，进而影响整个飞行器的气动性能。湍流输运系数变分项需要对湍流模型方程进行变分推导，先只考虑层流输运系数变分项的贡献，其中层流黏性系数μlam和层流热传导系数Klam的变分项可以通过Sutherland公式和层流Prandtl数Prlam计算：

(10)

式中：下标lam表示层流参数；T为温度；μref和Tref为Sutherland公式中的参考黏性系数和参考温度；Sref为Sutherland公式中的常数；cp为定压比热。将式(8)代入式(5)的第3项，采用Gauss积分定理和壁面无滑移条件简化可以得到：

(11)

式中：u为速度矢量；φ为与动量对应的伴随变量；S∞为远场边界；K为热传导系数；Σφ和η为中间变量；ψρE为与能量对应的伴随变量；μ为黏性系数;∂n表示法向梯度，式(11)中最后一项体积分项即为与NS方程对应的伴随方程，设ψρ为与密度ρ对应的伴随变量，Lagrange乘子Ψ=[ψρ, φ, ψρE]是与流场守恒变量相对应伴随变量。其中：

(12)

式中：h为焓值。推导式(11)的目的是建立目标函数与壁面变分的导数关系，进一步推导需要对NS方程的壁面边界条件进行变分运算，建立目标函数变分δJ′和壁面变分δSw的关系。

1.2 远场边界条件

对于气动力目标函数，在远场边界可以忽略黏性的影响，式(11)第2项表示在远场边界Γ∞上的积分。高超声速计算通常采用小远场边界，伴随方程边界条件需要结合流动控制方程边界条件进行分析。伴随变量在远场边界需满足

ΨT(Ac·n)δU=0

(13)

将式(13)展开，写为初始变量的形可得

(14)

式中：un为边界法向速度分量；γ为比热比。

对于超声速入口边界，边界流场变量密度、速度和压力均由来流计算条件指定，此时

(15)

已满足式(13)，不需要对伴随变量施加边界条件。

对于超声速出口边界，边界流场变量均由流场内点决定，此时

(16)

伴随变量需满足

(17)

在超声速出口边界Jacobian矩阵均为非零特征值，此时伴随方程边界条件为

(18)

对于亚声速出口边界，需要给定出口背压，此时：

(19)

伴随变量需满足

(20)

求解得

(21)

式中：c为声速；ψρE可由流场内点插值决定。

1.3 壁面边界条件

等温壁面上流动控制方程的边界条件为

(22)

式中：Tw为给定的壁面温度。线化处理后的变分形式为

(23)

将式(23)代入式(11)，为消除包含流场变量的积分项，得到连续伴随方程在等温壁面的边界条件：

(24)

式(11)可以转化为

(25)

式(25)中积分项表示气动力目标函数关于局部壁面变分的灵敏度：

giso=μn·Σφ·∂nu+K∂nψρE∂nT

(26)

式中：下标iso表示等温壁面。

对于绝热壁面条件，流动控制方程的边界条件为

(27)

线化处理后的变分形式为

(28)

式中：表示壁面切向梯度。将式(28)代入式(11)，得到伴随方程在绝热壁面的边界条件：

(29)

气动力目标函数在绝热壁面关于局部壁面变分的灵敏度为

(30)

式中：下标adi表示绝热壁面。

1.4 数值方法

采用基于非结构网格的有限体积法对流动控制方程和连续伴随方程进行求解。对于流动控制方程，对流项离散格式决定了流场波系结构的求解精度，对流项离散格式在强激波间断等位置需要具备强鲁棒性和低耗散性，在保证计算收敛的同时提高激波捕获精度。为此采用MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)格式[14]对界面变量进行重构，网格界面上对流项采用稳定性好、计算精度高的Roe格式[15-16]进行离散：

(31)

式中：下标ij表示网格控制体j对i的贡献，下标R、L表示重构后界面两侧流场变量对应的值；下标Roe表示界面两侧流场变量Roe平均后对应的值；PRoe为Roe平均后特征向量构成的矩阵；ΛRoe为Roe平均后对应的特征值。为在强激波和声速点位置满足熵条件，对Roe平均特征采用Harten-Yee型熵修正[17]：

(32)

式中：ε*为常数，通常取0～0.4；cRoe为Roe平均的声速。

对于连续伴随方程，为对流项，为黏性项，作为源项进行计算。伴随方程对流项与流动控制方程对流项具有共轭关系，离散时与流动控制方程面临相同的问题，伴随方程对流项在离散格式上与流动控制方程保持一致有利于提高伴随方程的求解精度[18]，伴随方程对流项的离散格式可以基于流动控制方程对流项的Jacobian矩阵进行构造：

(33)

式中：为伴随方程对流项；为对流项的第i个分量。式(33)采用一阶近似，同时应注意到伴随方程对流项在界面上不具有守恒性。流动控制方程和伴随方程黏性项采用中心差分格式离散，构造得到求解矩阵后，采用GMRES(Generalized Minimum RESiduals)方法[19-20]进行求解。

1.5 优化框架

单次优化循环的主要计算耗时集中在流动控制方程和连续伴随方程的求解，连续伴随方程求解计算量与流动控制方程相当，通过一次流动控制方程和连续伴随方程的求解可以得到目标函数关于所有设计变量的梯度，因此伴随方法具有计算量与设计变量数目无关的优势，适用于多设计变量优化问题。

图1给出了基于伴随方法的梯度优化流程，首先对优化对象建立几何参数化模型，采用FFD(Free Form Deformation)方法[21-22]对优化对象壁面进行参数控制。FFD方法可以较少的设计变量实现对变形物体的光滑连续控制；对于三维复杂外形可以方便地建立参数化模型；具备较强的局部变形控制能力[23]。通过网格变形方法将壁面变形传递到整个计算网格。优化对象气动性能通过求解RANS(Reynolds-Averaged Navier-Stokes)方程进行评估，伴随求解器读取流场变量对连续伴随方程进行求解得到目标函数壁面灵敏度，根据FFD方法将壁面灵敏度映射为目标函数关于设计变量的梯度。梯度优化算法采用SQP(Sequential Quadratic Programming)算法[24]，SQP算法根据目标函数梯度信息开展优化循环。

图1 优化流程

Fig.1 Optimization procedure

2 流场求解器验证

为验证流场求解器精度，对Sanger飞行器在高超声速来流条件下的气动性能计算[25]，并与AHL3D(Airbreathing Hypersonic Lab 3D)计算结果进行对比。AHL3D是中国空气动力研究与发展中心开发的大规模三维数值模拟软件，经过了大量研究和工程应用验证[26]。AHL3D计算采用结构网格如图2所示，结构网格分为49个块，第1层网格高度为10-6，网格量约为1 080万，对流项离散采用AUSM格式，湍流模型采用k-ω 剪切应力传输(Shear Stress Transter,SST)模型。求解器采用非结构混合网格，壁面边界层采用三棱柱网格，第1层网格高度为10-5，流场中间采用四面体网格，对流项离散为Roe格式，湍流模型采用带可压缩修正的SA模型，计算条件为来流马赫数6.0、高度35 km，壁面条件为1 000 K等温壁面。图3给出了不同攻角下升阻力系数对比，本文求解器得到的升阻力系数相对AHL3D偏小，阻力系数相差在5.0%以内，升力系数相差在3.0%以内。

图2 Sanger飞行器结构网格

Fig.2 Structured grids of Sanger vehicle

图3 Sanger飞行器升阻力系数对比

Fig.3 Comparison of lift and drag coefficients of Sanger vehicle

为对流场求解器网格无关性进行研究，通过调整非结构网格的表面网格尺寸和空间网格增长率，得到了如图4所示的3套网格，网格量分别为267万、494万和1 108万，来流攻角为4°。表1给出了3套网格升阻力系数计算结果，其中相对误差为相对网格3计算结果的误差，升力系数最大相对误差为1.9%，阻力系数最大相对误差为0.8%。

图4 Sanger飞行器非结构网格

Fig.4 Unstructured grids of Sanger vehicle

表1 Sanger飞行器气动力系数网格灵敏度

Table 1 Grid sensitivity of Sanger aerodynamic coefficients

3 伴随求解器验证

采用NACA64A-204翼型对伴随方法求解结果进行验证，计算网格采用结构网格的形式(图5)，采用FFD方法对翼型进行参数化控制(图6)，设计变量为FFD控制点沿法向的位移量，来流条件和计算设置与第2节中相同，攻角为4°。

图5 NACA64A-204翼型计算网格

Fig.5 Computational grids of NACA64A-204 airfoil

图6 NACA64A-204翼型FFD模型

Fig.6 FFD model of NACA64A-204 airfoil

目标函数选择翼型升阻比，图7给出了壁面灵敏度积分值和伴随变量残差的收敛历程，图中为灵敏度在整个壁面的积分值，ψρ u和ψρ v分别为与守恒变量ρu和ρv对应的伴随变量，其中v为初始变量中的速度分量,计算过程伴随变量残差一致收敛，经过3万次迭代求解，翼型壁面灵敏度积分值基本稳定，伴随变量残差降低了2个数量级以上。图8给出了伴随变量ψρE分布云图，由于伴随方程与流动控制方程的共轭关系，伴随变量的分布呈现与流场相反的特征，翼型前缘呈现出明显的间断分布，说明所采用的离散格式较好地捕捉到了伴随流场结构。图9给出了伴随方法得到的升阻比关于设计变量梯度与有限差分结果的对比，有限差分步长分别取为1×10-3、5×10-4、1×10-4，不同步长得到的有限差分结果一致性较好，同时前向差分(1×10-3)和后向差分(-1×10-3)结果符合较好。有限差分法和伴随方法在第0、11个设计变量处出现了明显的差异，在其余设计变量处符合较好，这是因为第0、11个设计变量对应的是翼型前缘，翼型前缘头部激波强度较大，伴随方程及其对应的边界条件在推导过程中采用了一阶线性近似，激波位置剧烈变化的流场参数会对计算精度产生一定影响。

图7 NACA64A-204翼型伴随方程收敛历程

Fig.7 Adjoint equation converence histories of NACA64A-204 airfoil

图8 NACA64A-204翼型ψρE云图

Fig.8 NACA64A-204 airfoil ψρE contour

图9 NACA64A-204翼型伴随方法与有限差分法梯度对比

Fig.9 Gradient comparison of adjoint and finite difference of NACA64A-204 airfoil

4 Sanger飞行器机翼优化

采用伴随优化框架对Sanger飞行器进行验证和优化研究，Sanger飞行器为德国二级入轨空天运输系统的一级飞行器，采用后掠三角翼兼顾中低速气动性能。机翼为影响整个飞行器气动性能的主要部件，采用FFD方法对飞行器机翼进行参数化控制，上下表面共设置64个设计变量(图10)，设计变量取值范围以1～8设计变量所在剖面翼型厚度为基准，为避免机翼前缘过度变薄，控制机翼前缘变形的设计变量的取值范围为基准的20%，其余设计变量的取值范围设定为基准的50%。计算网格采用第2节中的网格2(图4(b))。

图10 Sanger飞行器机翼FFD控制模型

Fig.10 FFD control model of Sanger vehicle wing

为对连续伴随方法进行应用验证，选择Sanger飞行器巡航状态进行优化，巡航来流状态为马赫数6.0、高度35 km、攻角4°，计算设置与第2节中相同，目标函数选择升阻比。图11给出了壁面灵敏度积分值和伴随变量残差的收敛历程，图中ψρw为与守恒变量ρw对应的伴随变量，计算过程伴随变量残差一致收敛，经过2万次迭代求解，翼型壁面灵敏度积分值基本稳定，伴随变量残差降低了2个数量级以上。图12给出了伴随变量ψρE在对称面的分布云图，可见采用的数值离散格式较好地捕捉到了伴随流场的波系结构。图13 给出了升阻比壁面灵敏度分布，红色区域表示壁面内凹时飞行器升阻比增大，蓝色区域反之。图14给出了伴随方法计算得到的升阻比关于设计变量梯度与有限差分结果的对比，由于设计变量数目较多，只选择同一控制面上的设计变量9～16(机翼下表面)、41～48(机翼上表面)进行有限差分验证，差分步长选择与NACA64A-204算例相同，可见不同步长得到的有限差分结果一致性较好，同时前向差分和后向差分结果符合较好。从对比结果看，相对二维情况下此时伴随方法与有限差分法结果在数值上的差异更明显，这是由于三维情况下流场激波结构更加复杂，同时机翼下表面气流受到的压缩更强，这种差异在设计变量9～16上更加明显；但两种方法得到的结果变化规律和正负分布一致，能够为梯度优化算法提供正确的优化方向。

图11 Sanger飞行器伴随方程收敛历程

Fig.11 Adjoint equation converence histories of Sanger vehicle

图12 Sanger飞行器对称面ψρE云图

Fig.12 Symmetric ψρE contour of Sanger vehicle

图13 升阻比壁面灵敏度

Fig.13 Lift-drag ratio surface sensitivity

图14 Sanger飞行器目标函数梯度对比

Fig.14 Comparison of objective function gradients of Sanger vehicle

图15给出了飞行器升阻比和3个站位翼型最大厚度的变化历程，优化循环共迭代29步，3个站位的位置分别为0.57(站位1)、0.72(站位2)和0.86(站位3)展宽位置，优化过程中站位1处翼型最大厚度变化较小，站位2和站位3处翼型最大厚度出现明显下降。飞行器升阻比在第11～20迭代步出现震荡收敛，这是由于优化历程趋近收敛时，伴随梯度精度的不足在一定程度上对优化历程产生影响，此时SQP算法通过减小试探步长在一定程度上克服了这种影响，从而导致飞行器升阻比出现上下震荡。表2给出了优化前后飞行器升阻比、升阻力系数和3个站位翼型相对厚度的变化，优化后飞行器升阻比增加5.0%，升力系数增加1.8%，阻力系数减小2.9%，3个站位翼型厚度分别减小1.8%、23.6%和25.4%。

表2 优化结果

Table 2 Optimization results

图15 优化历程

Fig.15 Optimization histories

图16给出了优化后飞行器机翼构型变化，优化后机翼上表面局部出现了小幅隆起，机翼前缘变薄，下表面出现明显的内凹，机翼构型变化与图13 中壁面灵敏度的分布一致，验证了FFD方法具备较强的局部变形控制能力。图17给出了3个站位剖面翼型的变化，机翼下表面内凹使得剖面翼型前缘厚度减小、弯度增加。

图16 优化后Sanger飞行器机翼构型

Fig.16 Optimized configuration of Sanger vehicle wing

图17 3个站位机翼构型变化

Fig.17 Configuration changes of vehicle wing at 3 stations

图18和图19分别给出了优化前后3个站位压力系数沿流场截面和翼型表面的变化，优化后机翼前缘迎风面积减小，头部激波强度减弱，减小了飞行器阻力；气流经过头部激波后出现明显膨胀，在内凹拐点处经过二次激波压缩，大幅增加了机翼下表面中后部压力，显著提高了飞行器的升力；优化后站位2上表面压力由于机翼隆起出现了小幅度下降，但机翼上表面压力变化较小，下表面压力变化是飞行器升阻比提升的主要因素。图20给出了优化前后机翼下表面的流场变化，从三维流动分析来看，优化后机翼前缘压力明显减小，紧贴机翼下表面的二次激波封闭了高压气流，减少了机翼前缘高压气流的泄流，机翼下表面高压区范围增大，达到了减阻增升的效果。

图18 3个站位截面压力系数分布

Fig.18 Pressure coefficient distribution of 3 station sections

图19 3个站位压力系数分布变化

Fig.19 Pressure coefficient distribution changes of 3 stations

图20 机翼附近流场变化

Fig.20 Changes of flow field near wing

5 结 论

本文构建了基于连续伴随方法的高超声速飞行器高精度气动优化方法，通过对二维翼型和Sanger飞行器机翼优化验证和应用得到以下结论。

1)基于熵修正的Roe格式构造了伴随方程对流项离散格式，在二维和三维伴随流场求解过程中伴随变量残差一致收敛，较好地捕捉到了与初始流场共轭的伴随流场波系结构，表明采用的数值格式具有较好的鲁棒性和低耗散性。

2)基于构造的数值格式，建立了适用于高超声速流动的连续伴随求解器，通过二维翼型算例和Sanger飞行器机翼优化算例的验证分析可知，在高超声速来流条件下构造的连续伴随求解器能够较好地给出气动力目标函数的梯度信息，但计算精度一定程度上仍受到强激波的影响。

3)FFD方法具备较强的局部变形控制能力，与基于连续伴随方法的梯度优化相结合，实现了针对高超声速飞行器的高精度优化。

4)Sanger飞行器机翼经过优化后，机翼前缘内凹，弯度增加，前缘迎风面积较小，产生的二次激波压缩使机翼下表面压力显著增加，飞行器升阻比增加5.0%，达到了减阻增升的效果。

5)经过本文的验证和应用研究可知，连续伴随方法是实现高超声速飞行器高精度气动优化的可行方法。