0 引言

以超燃冲压发动机为动力的吸气式高超声速飞行器是未来实现空天飞行和全球快速到达的重要途径。进气道作为压缩部件，其工作效率是影响冲压发动机工作状态和性能的重要因素。研究表明，对于Ma=5~7的碳氢燃料发动机，进气道压缩效率每提高1%，发动机比冲可提高3%~5%^[1]。高超声速进气道是机体和发动机强耦合部件之一，已经试飞成功的美国X-43A和X-51A飞行器分别采用了二维曲面压缩和内转进气道。这类进气道是基于无黏特征线反设计得到，在真实飞行条件下，受黏性和激波边界层干扰的影响，进气道的性能恶化。对于二维简单构型进气道可以采用黏性修正方法提高进气道性能，但对于三维复杂构型进气道，黏性修正方法很难得到理想的结果。在结构设计上进气道构型还应满足吸气式高超声速飞行器一体化设计的几何约束^[2]。从提高进气道性能和一体化设计的角度考虑，进气道优化设计是吸气式高超声速飞行器技术研究的重要内容。

以遗传算法为代表的随机类算法是目前高超声速飞行器优化设计所采用的主要方法。遗传算法鲁棒性强，可以得到目标函数的全局最优解，但遗传算法在优化循环中搜索面广，需要评价每个个体的优劣，计算开销较大，对于多设计变量的优化问题存在困难。为了提高优化效率，在工程应用中通常采用简化分析模型或通过少量的样本点构建代理模型，但针对多设计变量问题构建精确的代理模型仍需要较大的计算开销。

另一类优化方法是基于目标函数梯度的局部寻优算法，如拟牛顿方法和逐步二次规划方法，这一类算法是对目标函数采用局部二次多项式近似，根据目标函数的梯度信息确定寻优方向，逐步向目标函数的局部最优逼近。梯度信息的获取可以通过有限差分法，有限差分法需要分别计算目标函数对每个设计变量的梯度，计算开销随设计变量的数目增长。另一种方法是伴随方法，伴随方法通过构建流动控制方程的伴随方程，通过一次伴随方程的求解可以获得目标函数对所有设计变量的梯度值，计算开销与设计变量数目无关，伴随方程的求解计算量与流动控制方程大致相当。

伴随方程最早由Pironneau在控制理论中引入^[3]，Jameson从应用的角度出发，推导了势流方程和Euler方程的伴随方程，首次给出了基于伴随方法的二维翼型反设计算例^[4]，随后Jameson^[5-7]、Anderson^[8]和Giles^[9-10]等对连续伴随方法和离散伴随方法的推广以及在非结构网格上的应用都做出了大量的工作，目前伴随方法在CFD误差评估^[11]、不确定度分析^[12]和自适应网格^[13]等问题中都有重要应用。

在气动优化研究中，伴随方法主要应用于翼型、叶轮机械的优化^[14]。目前，国外研究人员已经开展了伴随方法在高超声速优化问题中的应用探索，Copeland给出了考虑化学非平衡NS方程的连续伴随方程的推导和求解，并针对再入航天器热流优化问题对比了伴随方法和有限差分法计算得到的梯度值^[15]；Kline采用连续伴随方法针对单楔压缩的二维和准三维高超声速进气道优化问题进行了研究^[16-17]。

伴随方法在高超声速进气道优化中的应用处于初步尝试和验证阶段，为了进一步探索伴随方法的应用潜力，本文选取二维曲面压缩无黏进气道构型，采用FFD方法对进气道外压缩面实现参数化控制，在考虑黏性的情况下以进气道出口流量为目标函数求解连续伴随方程，采用三套网格验证网格密度对外压缩面壁面灵敏度的影响，通过与有限差分法对比验证连续伴随方法计算结果的可靠性。选取外压缩面为优化对象，采用BFGS算法优化进气道出口流量，验证伴随方法在高超声速进气道优化中的适用性。

1 二维曲面压缩进气道

二维曲面压缩进气道设计采用特征线理论^[18]，根据求解问题类型的不同，进气道分为4个区进行反设计：前缘激波依赖区(A区)、主压缩区(B区)、末端激波依赖区(C区)和稳定区(D区)(图 1)。A区的求解依赖于来流条件、初始压缩角和壁面压升规律，B区的求解依赖于A区出口参数分布和壁面压升规律，C区的求解依赖于B区出口参数分布并在肩点c处实现反射激波消波，D区对气流方向进一步梳理。设计时，给定的设计马赫数为6.3，设计高度为27 km，初始压缩角为6°，壁面压升规律选择双曲正切函数。设计得到的进气道总长为L=849.8 mm，捕获高度为H=139.1 mm。图 2给出了特征线方法和Euler方程计算得到的进气道流场马赫数分布，Euler方程计算得到的流场结构与特征线方法符合较好，在Euler方程计算结果中，内压缩段肩点c处出现了微弱的反射波。

2 伴随优化方法 2.1 伴随方程

连续伴随方法是基于流动控制方程推导而来，以进气道出口流量为目标函数，其可表示为进气道出口的积分:

$ J = \int_{{\mathit{\Gamma }_e}} g {\rm{d}}s = \int_{{\mathit{\Gamma }_e}} \rho \mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{n}}{\rm{d}}s $

(1)

定常流动控制方程的差分形式表示为:

$ \left\{\begin{array}{l} \boldsymbol{R}(\boldsymbol{U})=\nabla \cdot \boldsymbol{F}=\nabla \cdot \boldsymbol{A} \boldsymbol{U}=0 \quad \text { in } {\mathit{\Omega}} \\ \boldsymbol{v}=0 \quad \text { on } \quad S \\ \boldsymbol{v}=\boldsymbol{v}_{\infty} \quad \text { on } \quad {\mathit{\Gamma}} _{\infty} \end{array}\right. $

(2)

其中 R 为残差，U 为流场守恒变量，A 为Jacobi矩阵。优化问题可以等价为以流动控制方程为约束，求目标函数的极值。通过Lagrange乘子将流动控制方程引入目标函数:

$ {\mathscr{T}}=J-\int_{{\mathit{\Omega}}} \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{R}(\boldsymbol{U}) \mathrm{d} {\mathit{\Omega}} $

(3)

目标函数关于壁面的灵敏度$\partial \mathscr{T} / \partial S$表示了壁面法向变形对目标函数的影响量，可以用目标函数的变分δ${\mathscr{T}}$/δS进行近似:

$ \delta \mathscr{T}=\delta J-\int_{\mathit{\Omega}} \boldsymbol{\psi}^{\mathrm{T}} \delta \boldsymbol{R}(\boldsymbol{U}) \mathrm{d} {\mathit{\Omega}} $

(4)

运用变分原理将流动控制方程代入公式(4)可得:

$ \begin{aligned} \delta \mathscr{T}=& \int_{{\mathit{\Gamma}}_{e}} \frac{\partial g}{\partial \boldsymbol{U}} \delta \boldsymbol{U} \mathrm{d} s-\int_{{\mathit{\Gamma}}} \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n} \delta \boldsymbol{U} \mathrm{d} s-\\ & \int_{S} \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n} \delta \boldsymbol{U} \mathrm{d} s-\int_{S} \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n} \delta S \mathrm{d} s+\\ & \int_{{\mathit{\Omega}}} \nabla \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \delta \boldsymbol{U} \mathrm{d} {\mathit{\Omega}} \end{aligned} $

(5)

为了计算δ$\mathscr{T}$/δS，可以通过构建伴随方程及相应的边界条件来消除公式(5)中包含δU 的项:

$ \left\{\begin{array}{l} \nabla \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A}=0 \quad \text { in } \quad {\mathit{\Omega}} \\ \frac{\partial g}{\partial U}-\boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n}=0 \quad \text { on } \quad {\mathit{\Gamma }_\mathit{e}} \\ \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n}=0 \quad \text { on } \quad \Gamma \neq {\mathit{\Gamma }_\mathit{e}} \\ \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n}=0 \quad \text { on } \quad S \end{array}\right. $

(6)

上式即为流动控制方程(2)的伴随方程，ψ 称为伴随变量。通过构建伴随方程和相应的边界条件，目标函数的灵敏度就可以表示为:

$ \delta \mathscr{T}=\int_{S} \boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{A} \cdot \boldsymbol{n} \delta \boldsymbol{S} \mathrm{d} s $

(7)

在应用伴随方法时，首先对流动控制方程进行求解得到的流场变量，代入伴随方程，迭代求解得到伴随变量，将伴随变量代入公式(7)即得到进气道流量关于壁面的灵敏度信息。

2.2 曲面参数化和网格变形

在气动外形优化问题中，需要通过正交、完备的基函数对优化型面的设计空间进行参数化表示，可以减少优化过程中的设计变量，同时在优化循环中保证型面的光滑和连续。采用FFD方法对进气道外压缩面进行参数化控制，FFD方法通过构建包裹变形曲面的控制体，利用基函数建立变形曲面与控制点的映射关系:

$ \boldsymbol{X}(u, v, w)=\sum\limits_{i, j, k=0}^{l, m , n} \boldsymbol{P}_{i, j, k} B_{i}^{l}(u) B_{j}^{m}(v) B_{k}^{n}(w) $

(8)

其中 X 为变形曲面上点的坐标，P 为控制点的坐标，B_i^l(u)=C_lⁱuⁱ(1－u)l^－i为l阶的Bernstein多项式，通过映射关系公式(8)利用控制点的移动来控制型面的变形。图 3给出了进气道外压缩面的FFD控制体的示意图，设置FFD控制体的维度为11×2，保证外压缩面在起始和终点位置具有连续的2阶导数，设计变量选择(0, 1)~(0, 9)9个控制点的y坐标值，其余控制点位置固定不变。从图上可以看出，在变形后外压缩面保持了光滑和连续。

网格变形采用弹性比拟法，将计算网格等效为弹性体，以外压缩面的位移为边界条件来求解网格点的变形量，通过求解线弹性Poisson方程得到各个网格点的位移量:

$ \left\{\begin{array}{l} \nabla \cdot \boldsymbol{\sigma}=\boldsymbol{f} \text { in } {\mathit{\Omega}}\\ \boldsymbol{\sigma}=\lambda {Tr}(\boldsymbol{\varepsilon}) \boldsymbol{I}+2 \mu \varepsilon \\ \boldsymbol{\varepsilon}=\frac{1}{2}\left(\nabla \boldsymbol{u}+\nabla \boldsymbol{u}^{\mathrm{T}}\right) \\ \boldsymbol{u}=\boldsymbol{u}_{\mathrm{S}} \quad \text { on } \mathrm{S} \end{array}\right. $

(9)

其中Tr(ε)为ε的迹，λ=νE(1-2ν)/(1+ν)，μ=2E/(1+ν)，E为弹性模型，ν为泊松比。在求解网格变形过程中，边界层和局部加密区域网格由于网格纵横比大、体积小，这些网格容易过渡扭曲产生负体积，在求解过程中采用随网格体积成反比变化的弹性模量，增强这些网格的刚度，保持原有的网格质量。

3 外压缩面灵敏度分析

进气道的初始构型为采用特征线法设计的无黏构型，在黏性情况下由于受边界层挤压的影响，唇口处会产生溢流，外压缩面形状影响到前体流场波系结构进而影响进气道的捕获流量。计算构型前缘和唇口进行倒圆(倒圆直径1 mm)，采用三套网格(如图 4所示)对目标函数进气道流量关于外压缩面灵敏度的影响因素进行分析。第一套网格(Grid 1)壁面网格第一层高度为0.1 mm，网格总数为3.9万；第二套网格(Grid 2)壁面第一层网格高度为0.01 mm，流场中间网格分布与第一套网格基本保持一致，网格总数为4.6万；第三套网格(Grid 3)在第二套网格的基础上在对流场中间网格进行了加密，壁面网格基本保持不变，网格总数为6.5万。

在计算外压缩面灵敏度时，伴随方程考虑了湍流黏性项的影响。流动控制方程和伴随方程求解过程中，无黏通量计算采用Roe格式，湍流模型采用SA模型。

图 5给出了在三套网格上流动控制方程(Flow)求解过程中进气道流量的收敛历程和伴随方程(Adjoint)求解过程中外压缩面总灵敏度的收敛历程。对于同一套网格，流动控制方程和伴随方程的收敛历程保持一致；三套网格得到的进气道出口流量相差很小，Grid 1计算得到的外压缩面总灵敏度相对偏小；三套网格迭代60 000步后流动控制方程和伴随方程都基本收敛。

图 6给出了三套网格计算得到的外压缩面灵敏度的分布，在进气道前缘处受前缘曲率的影响，灵敏度出现异常峰值。为了避免对优化过程产生误导，在伴随方程求解过程中设置光顺因子(Smooth Ratio，SR)，将前缘附近的灵敏度进行光顺处理。对比Grid 2得到的两条灵敏度曲线(Grid 2 SR=0.01和Grid 2 SR=0.03)，当光顺因子较小时(SR=0.01)，前缘附近出现了明显的异常峰值，当设置光顺因子SR=0.03，前缘峰值消失，同时正常峰值(x=0.47 m处)也有所降低。

各条灵敏度曲线在x=0.47m处均出现峰值，分析进气道流场，该处产生的压缩波与入射激波在唇口处交汇，直接影响了进气道的捕获流量。在峰值之后产生的压缩波入射进入进气道内部，对捕获流量不再产生影响，所以当x>0.47 m后灵敏度迅速下降到0。Grid 2和Grid 3计算得到的灵敏度分布基本一致，Grid 1计算得到峰值相对较小，说明壁面灵敏度结果对流场中间网格分布不敏感，对边界层网格分布依赖性较强。

为了对灵敏度的合理性进行验证，根据公式(11)将灵敏度分布映射到9个设计变量，得到目标函数关于设计变量的梯度。同时根据有限差分法，分别在每个设计变量上施加差分量，计算目标函数对该设计变量的梯度，差分步长分别取1×10^-3m，5×10^-4m和2.5×10^-4m，计算得到的结果如图 7所示。在不同步长下差分方法得到的结果相互吻合，说明有限差分结果具备收敛性。对于伴随方法，Grid 2和Grid 3得到的结果略高于差分方法，变化规律一致，第一套网格得到的结果明显偏小，说明伴随方法可以得到有效的灵敏度信息。

4 优化结果

根据以上建立的参数化几何模型和伴随方程求解方法，以进气道流量为目标，对进气道外压缩面进行优化，计算网格选择第二套网格，优化算法采用BFGS方法，BFGS方法是一种改进的拟牛顿法，是利用目标函数梯度求解无约束优化问题的最有效的一类方法。图 8给出了进气道流量的优化历程, 在第3次优化循环后流量达到稳定，流量由12.37 kg/s增加到13.15 kg/s，提高了6.3%。

图 9给出了优化前后的流场马赫数分布。在优化前，对于无黏进气道构型，在黏性情况下气流受到边界层排挤，唇口处产生了溢流。优化后，入射激波位置上移实现封口，消除了溢流，基本实现了全流量捕获。图 10给出了优化前后外压缩面型线和压力系数的分布变化，优化后外压缩面前部出现了略微的下凸，在中后部出现明显的内凹，与灵敏度的分布规律一致，最大法向变形为5.6 mm。优化后外压缩面压力梯度整体增大，外压缩面压缩能力提升。表 1给出了优化前后进气道出口流量系数φ，质量平均马赫数Ma_exit，压升比p_exit/p_∞和总压恢复系数σ的变化，优化后进气道出口马赫数提高约2.4%，由于外压缩面压缩能力的提升，进气道出口压升比提高约5.3%，外压缩面以等熵压缩为主，整体上减少了进气道的总压损失，总压恢复系数提高了约10.6%，进气道整体性能得到提高。

5 结论

本文针对连续伴随方法在二维高超声速进气道优化设计中的应用开展了研究。进气道构型采用基于特征线法设计的曲面压缩进气道。为了分析壁面灵敏度的影响因素，采用了三套疏密不同网格进行计算对比分析。对于目标函数灵敏度，采用有限差分结果进行对比验证。基于所采用的连续伴随方法和FFD参数化方法，以进气道外压缩面为优化对象，采用BFGS算法对进气道流量进行优化，得到了以下结论：

1) 以进气道流量为目标函数，伴随方法可以有效得到外压缩面灵敏度信息，灵敏度对边界层网格分布有较强的依赖性。

2) 采用BFGS算法对进气道流量优化，优化后唇口处实现激波封口，流量提高6.3%，外压缩面最大法向位移为5.6mm，外压缩面压缩能力增强，进气道整体性能得到提高。

3) 从优化结果来看，伴随方法可有效应用于高超声速进气道优化。