基于共轭方程法的跨声速机翼阻力优化设计

张野平， 段卓毅， 张彦军

（西安飞机设计研究所，陕西 西安710089）

0 引 言

在飞行器设计中，翼型、机翼、全机的气动设计技术历来被视为各国航空领域的核心技术，传统的设计方法主要有间接方法、余量修正方法几种类型。间接设计法中典型的有虚拟气体法和速度图法。余量修正方法中具有代表性的方法是Takannashi在19世纪80年代中期提出的一种余量修正方法［1］。这两类方法都计算量相对较小，但是自身存在不能直接控制气动特性或者解的适定性等诸多问题。

20世纪70年代初至今，随着计算机科学和数值方法的迅猛发展，计算流体力学（CFD）和数值优化技术取得了长足的进步，基于CFD的气动外形优化设计技术应运而生。根据目前的计算能力而言，以遗传算法为代表的全局寻优方法由于搜索面广、计算量大，还只能用在比较简单的模型上；能够应用于工程实际的只有基于梯度计算的局部寻优方法。在基于梯度的气动优化设计方法中，对于有N个设计变量的设计问题，往往进行一次优化设计就要进行N＋1次流场求解，为获得梯度所需进行的流场计算次数与设计变量个数直接相关，当设计变量众多时计算代价难以接受。

20世纪80年代末期，A.Jameson提出的基于控制理论（共轭方程）的气动优化设计技术［2－4］成功地解决了基于梯度的优化方法中梯度获得困难、计算量巨大这一难题。该方法以偏微分方程系统的控制理论为基础，把物体边界形状作为控制函数，流场控制方程作为约束条件，在目标函数中引入Lagrange因子，变条件约束问题为无约束问题，通过求解流场控制方程及其共轭方程进行梯度求解，其计算量只相当于两倍的流场计算量，与设计变量数目无关，从而将设计问题转变成控制问题。目前该设计理论的应用领域不断拓宽，对该方法的深入研究具有重要的意义和工程价值。

1 流动控制方程

控制方程使用NS方程。在三维笛卡尔坐标系（x1，x2，x3）中，有速度分量为（u1，u2，u3），守恒形式的NS方程形式为：

其中w是状态矢量，fi是无粘矢通量，fvi是粘性矢通量。

在用结构网格进行求解时，需要将上述NS方程转化为计算坐标系（ξ1，ξ2，ξ3）下，其守恒形式转化为：

其中F＝Sf，F＝Sf，S＝J，J为坐标变换iijjviijvjij的Jacobin值。Sij的物理意义是包含控制体的边界的外法向面矢量，J为控制体体积的倒数。

2 伴随方程

引入目标函数表示为：

各变量的含义如下：

w：流场参数

S：变换矩阵

B：物面边界

D：计算体积空间

M：物面边界变化对泛函的贡献

P：计算区域体积空间变化对泛函的贡献

很显然，就空间划分来说，目标函数受到的影响来自于物面边界的变化和计算空间单元体积变化两个部分。就物理上来说，泛函的值取决于流场参数w和变换矩阵S。

物面边界的变化最终将导致目标函数的变化。则可以物面边界作为控制函数，将流动控制方程作为约束条件在目标函数中引入，从而将约束问题转化为无约束问题，设计问题转化为控制问题。根据变分法则，设目标函数对于边界变形的变分为：

下标Ⅰ和Ⅱ分别表示由流动变量变化δw起的贡献和由坐标变换矩阵变化δS引起的贡献。为了将流动的控制方程作为约束条件引入到目标函数里面，从而将有约束问题转化成无约束问题，现在对定常情况下的流动控制方程进行变分处理：

在这里引入共轭变量Ψ，并对上式做全流场积分：＝0，假设“乘子”Ψ是可微分的，对上式做分部积分：

利用高斯定理，进行体积分到面积分的转化：

其中ni是面矢量。

将控制方程的积分形式的变分引入到目标函数的变分当中，用原来的变分减去现在的积分形式的控制方程的变分，再将δM、δP、δFi、δFvi带入，整理可得：

选择合适的共轭变量Ψ，使得δI和流动变量δw不再显性相关。目标函数的计算可直接通过坐标变换矩阵的变分得到，而不需要重新计算由扰动每个设计变量产生的流动变量变分δw。这样就可以在全流场区域D的积分当中首先得到全流场粘性伴随方程：

然后在边界积分区域B当中得到粘性伴随方程的物面边界条件：

可得

对于本文涉及的正向设计问题，目标函数一般表达为：

式中CD为压阻系数，CL为设计升力系数。

3 梯度求解

下面通过简单的变分计算问题说明改进过的Sobolev梯度求解方法［5］，使用该方法能够生成一系列连续的光滑外形，避免优化所得的外形连续性降阶。

选择y（x）使得I（y，y′）dx最小，y（x）在两端点的值分别为y（a）、y（b）。

使y（x）产生一个小的变量δy（x），可得：

定义梯度形式为：

内积形式为

可以看出

假设：δy＝－λg，λ＞0

由此可得：

当且仅当g＝0时，等号成立，即取最小值。注意到g是y、y′、y″的函数，即g＝g（y，y′，y″）。

以著名的Brachistrone（最速降线）问题为例：

并且：

每一步yn＋1＝yn－λngn都减小y的两阶光滑度。因此计算会导致越来越不光顺，将导致计算的不稳定性。

为了避免上述不利情况的发生，引入一种改进过的Sobolev内积：

这里ε是控制导数项权重的参数。如果定义一种梯度能够使得：

可以得到：

其中

则有每一步yn＋1＝yn－λnn。

这样就给出了一种改进以后δI的形式：

这种情况下yn＋1和yn具有一致的连续性，可以得到稳定的解。

使用这种求解梯度的方法，就可以只用网格点对模型进行参数化，而不必使用其它复杂的模型参数化方法。一般来说，使用每一个网格点来对模型进行参数化，即把每一个网格点当作设计变量，那么在三维机翼的情况下，设计变量将多达数千个，对于这么多的设计变量将产生一些问题，使用这样的设计空间，会导致所谓的病态问题。通过使用上述梯度求解方法，可以光滑控制表面形状，减弱了设计中的高频震荡，使病态问题得到了改善。

4 数值方法

计算使用剪切抛物变换得到的C－H型网格；湍流模型采用BL模型；在流动控制方程和共轭方程的数值求解中，采用格心格式的有限体积法进行空间离散，时间推进采用五步Runge－Kutta方法。

采用了目前常用的当地时间步长、隐式残值光顺和多重网格法加速计算收敛，同时加入人工粘性抑制数值震荡。

生成计算网格后直接使用网格点进行几何的参数化描述，每个网格点视作一个设计变量，通过文中所述的梯度求解方法可以避免优化所得外形的光滑性降阶；每完成一次优化计算，重新生成网格进行计算。

在求得梯度之后，使用最速下降法进行优化。

为保证设计所得的机翼在较大的范围内都具有良好的气动特性，将优化方法和Pareto阵面的概念结合，进行多目标气动优化［6］。

5 设计算例

优化设计构型为某飞机翼身组合体构型，仅对机翼部分进行优化设计。

为保证设计所得的机翼在较大的范围内都具有良好的气动特性，提出四个设计状态，如表1所示。

表1 设计状态Table 1 Design points

第一个设计点用以对爬升状态点进行控制；第二、三个设计点用以对机翼巡航状态进行减阻设计；第四个设计点用以对较高马赫数进行减阻优化，获得较高的阻力发散马赫数。

具体约束为机翼控制剖面的相对厚度、机翼的平面形状不改变，考虑到机翼的低速特性与机翼控制剖面的头部半径关系密切，另外加入对机翼控制剖面的头部半径的约束。

选取一副经过初步设计的机翼作为初始外形，采用大量实践验证的几何扭转角分布进行剖面扭转以保证良好的几何特性和翼尖失速特性，并且根据实际的剖面厚度要求对基本剖面形状进行厚度调整。初始机翼使用的控制剖面具有比较明显的超临界翼特征，这样配置所得的初始机翼具有较好的跨声速气动特性。

图1为翼身组合体表面网格和空间剖面网格示意图。

图2为优化前后机翼六个控制剖面外形对比。由于加入了剖面厚度、扭转角以及头部半径的约束，可以看出，各剖面厚度未出现明显变化，另外，优化后剖面头部半径均与初始机翼持平或者略大于初始机翼。

图3、图4为两个典型状态下优化前后剖面压力分布对比，图5、图6为多个状态下优化前后机翼表面压力云图，从中可以看出，较大马赫数下机翼上表面激波有明显减弱，这样可以大大降低波阻，提高阻力发散马赫数。图7为机翼优化前后CL＝0.55下阻力随马赫数变化曲线，图8为CL＝0.55下优化前后巡航效率因子随马赫数变化曲线。可以看出初始机翼阻力发散马赫数在0.76～0.77之间，经过15次优化后，阻力发散马赫数提升至0.78～0.79之间。巡航效率因子也较初始机翼有明显提高。

图1 表面网格及机翼展向剖面网格Fig.1 Surface mesh and spanwise section mesh

图2 机翼各展向位置剖面几何外形变化Fig.2 Variation of section shape at different spanwise positions

图3 Ma＝0.77，CL＝0.53剖面压力分布对比Fig.3 Pressure comparison of section，Ma＝0.77，CL＝0.53

图4 Ma＝0.79，CL＝0.5剖面压力分布对比Fig.4 Pressure comparison of section，Ma＝0.79，CL＝0.5

图5 Ma＝0.77，CL＝0.53机翼表面压力云图对比Fig.5 Comparison of pressure contours of wing，Ma＝0.77，CL＝0.53

图6 Ma＝0.79，CL＝0.5机翼表面压力云图对比Fig.6 Comparison of pressure contours of wing，Ma＝0.79，CL＝0.5

图7 CL＝0.55下阻力系数随马赫数变化曲线Fig.7 Variation of drag coefficient with Mach Number，CL＝0.55

图8 CL＝0.55下巡航效率因子随马赫数变化曲线Fig.8 Variation of cruise efficient with Mach Number，CL＝0.55

表2为四个设计点优化前后机翼主要气动力指标对比，可以看出各个设计点的阻力均有所减小，升阻比相应提高，大马赫数下优化效果更为明显，由于优化增大了机翼的后加载，因此机翼的俯仰力矩系数有相应的增加，四个设计状态下平均俯仰力矩系数增加0.02左右，在巡航状态Ma＝0.75～0.77附近，俯仰力矩系数约为CM＝－0.145附近，量值对于超临界机翼来说属于较为适中的范围，工程上可以接受。

表2 优化前后气动力对比Table 2 Comparison of aerodynamic performance before and after optimization

6 结 论

本文通过求解三维NS方程及其共轭方程对跨声速机翼进行气动优化，文中针对某三维机翼气动外形优化设计效果明显，算例证明本文发展的基于粘性NS方程共轭方法的跨声速机翼优化设计具有良好的优化效果，相比其它优化方法，该方法具有优化效率高、不受设计变量数目限制的明显优势，优化结果满足设计要求，具有很高的工程实用价值。

［1］TAKANASHI S.An iterative procedure for three dimensional transonic wing design by the integral equation method［R］.AIAA Paper 84－2155，1985.

［2］JAMESON A.Aerodynamic design via control theory［R］.AGARD－CP－463，1989.

［3］JAMESON A，REUTHER J.Control theory based airfoil design using the euler equations［R］.AIAA Paper 94－4272－CP，1994.

［4］JAMESON A，PRINCE N A，MARTINELLI L.Optimum aerodynamic design using the navier－stokes equations［R］.AIAA Paper 97－0101，1997.

［5］JAMESON A，VASSBERG J C，SRIRAM SHANKARAN.Aerodynamic－structural design studies of lowsweep transonic wings［R］.AIAA Paper 2008－145，2008.

［6］TANG Z L.Pareto front capture using deterministic optimization methods in multi－criterion aerodynamic design［J］.TransactionofNanjingUniversityofAeronautics＆Astronautics，2006，23（2）：81－86.