DOI: 10.11817/j.issn.1672-7207.2016.02.012

基于伴随方法的叶片三维气动外形优化设计

刘浩^{1, 2}，张雷¹，李霄琳¹

(1. 吉林大学 机械科学与工程学院，吉林 长春，130025；

2. 富奥汽车零部件股份有限公司 研发中心，吉林 长春，130013)

摘要：在计算流体动力学的基础上，通过引入离散黏性伴随方法和基于径向基函数的网格变形技术，构建一种高效的叶片气动优化设计系统。利用此优化设计系统对跨音速转子Rotor37的叶片进行优化设计。研究结果表明：在近最大等熵效率工况点下，优化后叶片流道激波强度和流动损失有效降低，转子的总压比提高0.352%，等熵效率提高0.339%，质量流量提高0.362%。与传统的叶片优化设计方法相比，此优化方法能够显示分析叶片外形与目标函数之间的关系，并有效地降低优化设计时间。

关键词：计算流体动力学；伴随方法；径向基函数；网格变形；优化设计

中图分类号：V231.3             文献标志码：A         文章编号：1672-7207(2016)02-0436-07

Three dimensional aerodynamic shape optimum design for blade based on adjoint method

LIU Hao^{1, 2}, ZHANG Lei¹, LI Xiaolin¹

(1. School of Mechanical Science and Engineering, Jilin University, Changchun 130025, China;

2. Research and Development Center, Fawer Automotive Parts Limited Company, Changchun 130013, China)

Abstract: Based on the computational fluid dynamics, a high-efficiency aerodynamic optimum design system of blade was constructed by using discrete viscous adjoint equation method and combining radial basis function mesh deformation technology. This optimum design system was applied to optimize the profile of an transonic rotor Rotor37. The results show that the shock wave intensity and flow loss decrease in the flow passage, and the total pressure ratio increases by 0.352%, the adiabatic efficiency by 0.339%, and the mass flow by 0.362% under the condition of near peak adiabatic efficiency. Compared with the traditional blade optimum design method, the proposed optimization system can analyse the relationship between blade shape and the cost function explicitly, and reduce the time of optimized design effectively.

Key words: computational fluid dynamics; adjoint method; radial basis function; mesh deformation; optimum design

对定转子叶片形状的进行优化是提升发动机性能的有效途径之一。随着计算机及计算流体动力学技术的发展，数值优化成为发动机叶片设计过程中必不可少的环节。目前，叶片气动优化设计方法主要有2种：基于搜索算法的全局优化算法和梯度基优化算法。全局优化算法主要有遗传算法、人工神经网络法和模拟退火算法等，这类方法需要较多的设计变量，并且要在整个设计空间内进行迭代寻优，因此，需要耗费巨大的计算资源，设计周期长。梯度基优化方法比全局优化算法计算量小，它是以目标函数对设计变量的梯度为寻优方向的搜索优化算法，高效地求解梯度是此优化方法的关键。传统的梯度优化算法采用有限差分法对梯度进行求解，当设计变量较多时，计算量仍然很大。为了提高梯度计算的效率，JAMESON^[1]以偏微分方程的控制理论为基础，提出了伴随优化设计方法。自NADARAJAH等^[2-3]成功地将伴随方法引入到航空领域，对翼型、机翼及机身进行气动优化设计至今，伴随方法已成为空气动力学中的研究热点，并在外流场的气动优化中得到了广泛的应用。近年来，在透平机械内流领域，伴随方法也逐渐得到关注，PAPADIMITRIOU等^[4-5]采用伴随方法研究透平机械中总压损失的问题并对其进行优化设计，WANG等^[6-7]将伴随方法用于多级透平机械叶片的优化设计中，张朝磊等^[8-9]也成功将其应用于叶栅的气动优化设计和反设计，由此伴随方法在透平机械中应用的有效性得到了验证。本文作者将离散伴随方法应用到发动机转子叶片的三维气动优化设计中，并结合径向基函数网格变形技术，建立了高效的叶片三维气动优化设计体系。首先介绍了伴随方法及N-S方程的伴随方程推导，其次引入网格变形方法建立优化体系，最后利用此优化体系对NASA Rotor37叶片进行气动优化设计，并分析了优化后叶片内部流场的变化。

1  伴随方法

伴随方法是基于偏微分方程系统控制理论的气动设计方法，通过求解伴随方程来获得流场目标函数对设计变量的敏感性关系。

1.1  伴随方法原理

在给定的流动条件下，目标函数L，流场变量ω和几何物理边界α，可以用以下函数表示：

               (1)

假设几何边界有微小的变化时，对应目标函数的变化量为

            (2)

流动控制方程R在流场区域D内表示为

               (3)

流场的变分为

           (4)

通过引入伴随变量，使目标函数对流场变量的变分系数为0，以此消除流场变分对目标函数梯度的影响。因此，引入拉格朗日因子ψ，则式(2)可表示为

     (5)

令

             (6)

则式(5)可简化为

         (7)

方程(6)称为伴随方程，方程(7)与流场变量ω无关，只与几何变化相关。伴随优化设计方法的优点是只要分别求解一次流动方程和伴随方程就可获得全部梯度信息。因此，当采用较多的设计变量时，目标函数梯度G的计算并不需要额外的流场计算，比传统的采用有限差分法计算梯度的方法要节约很多计算成本。

1.2  N-S方程的伴随方程推导

三维N-S方程在直角坐标系下，可表示为

           (8)

变换到计算域坐标系下，可表示为

           (9)

式中：为计算域坐标；f_i为无黏通量；f_vi为黏性通量；，，，，。

对于定常计算，N-S约束方程的变分为

           (10)

将其乘以拉格朗日因子ψ，并积分：

         (11)

利用分部积分，有

 (12)

目标函数与流场边界条件和状态有关，可以用积分形式将其表示为

      (13)

其变分为

    (14)

其中：

          (15)

           (16)

将目标函数的变分式(14)减去式(12)，则有

       (17)

消除式(17)与流场变分相关的量，可以得到只与几何边界相关的目标函数梯度。

可以得到N-S方程的伴随方程为

     (18)

对应的几何边界伴随方程为

       (19)

2  径向基函数网格变形

RENDALL等^[10]提出了一种基于径向基函数的网格变形方法，并在此基础上拓展了计算效率更高的缩减控制点的径向基函数网格变形技术^[11-12]。径向基函数是一种灵活的多维空间插值方法，与B样条和多项式插值方法相比，它可以实现对非结构型网格的插值。

径向基函数是一个取值与空间距离相关的函数，在径向基函数网格变形中的基函数以欧式距离为变量，可以分为全局型基函数(global basis function)和紧支型基函数(compact basis function)，全局型基函数的值随着与中心点的距离增加而增大，所有函数值均不为0；紧支型基函数，其函数值随着与中心点的距离增大而减小，当距离达到紧支半径时，函数值为0。网格变形时，要使中心点的移动只影响其周围一定范围的网格，因此RBF网格变形一般使用紧支型基函数，表1列出了常用的紧支型基函数，其中变量ξ为欧式距离与紧支半径的比值，当ξ＞1时，函数值为0，即中心点的移动不会影响紧支半径以外的点。

         (20)

用于网格变形的插值函数是由一个径向基函数和多项式h组成，表达式如下：

表1  常用紧支型基函数

Table 1  Common compact basis function

        (21)

式中：s(x)为插值函数；为径向基函数的一般形式；x_i为第i个中心点的位置；x为空间任意一点的位置矢量；为空间任意一点到第i个中心点的距离；与第i个中心点相对应的权重系数。

对于网格变形，需要x，y和z这3个方向的插值函数。选择系数和多项式h，使插值函数在中心点处，满足

        (22)

同时对所有一阶多项式q，应满足

            (23)

多项式h的阶数由所选的径向基函数决定，若选择的径向基函数阶数n≤2，则多项式h(x)可表示为

       (24)

径向基函数的系数，线性多项式h(x)的系数β，可通过求解下式获得

         (25)

其中：M为插值矩阵。

   (26)

P为约束矩阵用于平衡多项式对插值函数的影响，由中心点坐标和单位向量组成：

         (27)

则

      (28)

           (29)

其中：

          (30)

通过上述公式可以确定空间节点在各个方向上的位移v_x，v_y和v_z：

          (31)

         (32)

          (33)

基于径向基函数网格变形的重要特点是只需要单独的移动网格节点，分别插值各个方向上的位移，而不需要考虑网格之间的连续性，这个特性使得径向基函数适用于结构型和非结构型网格的变形。

3  优化设计

本文的研究对象为NASA的跨音速转子Rotor37，流场计算采用有限体积法求解三维定常RANS方程，湍流模型采用Spalart-Allmaras模型，对流项采用二阶精度的迎风格式进行离散。网格划分为HOH结构型网格，网格数量约为220万，采用较多的网格是为了保证能够在流场计算中获得与目标函数相关的流场梯度信息。

3.1  优化算法

在工程设计最优化问题中最常见的是搜索算法，梯度下降法也称最速下降法^[13]是求解无约束优化问题最简单的方法之一，是一种局部搜索寻优算法。

最速下降法是沿梯度下降的方向搜索极值，算法流程如下：

1) 给定初始迭代点x₀，确定终止准则：终止误差ε和最大迭代步数k_max。

2) 确定下降方向p_k：。

3) 确定搜索步长α_k：满足。

4) 新的迭代点：。

5) 确定是否满足或者达到最大迭代步数k_max，若达到收敛准则，最优点为x_k+1，则终止迭代；否则令k=k+1，转至步骤2)继续搜索。

3.2  优化流程

优化设计目标函数L为转子的总压比最大，设计变量d为基于RBF网格变形的中心点，寻优算法采用最速下降法。

优化设计流程如图1所示，基本步骤总结如下：

1) 根据几何模型，划分网格，面网格用x_surf表示，体网格用x表示。

图1  优化设计流程图

Fig. 1  Flow chart of optimization design

2) 选择流场计算模型，设置边界条件，进行流场计算，得到流场解w。

3) 基于流场计算结果，利用伴随方法计算目标函数对叶片表面网格的梯度。

4) 根据目标函数对叶片表面的梯度信息，选择网格变形中心点的位置和数量。

5) 计算目标函数对网格变形中心点的梯度。

6) 根据中心点梯度和坐标，通过最速下降法，计算得到新的中心点坐标，以驱动流体网格变形。

7) 对变形后的网格，重新进行流场计算，在步骤5)和6)之间循环，直到优化目标函数残差收敛。

上述优化过程的关键是对梯度的求解，文中通过以下公式进行简要说明。

设计变量d与面网格x_surf的关系可表示为

              (34)

则体网格x为

              (35)

流场解w为

             (36)

根据式(1)目标函数L可表示为

      (37)

由式(7)可得梯度：

         (38)

目标函数对网格变形中心点的梯度可由下式获得：

       (39)

图2所示为目标函数对叶片外形的梯度分布云图，它反映了目标函数对叶片面型变化的梯度信息，可以准确锁定优化区域。由图2可知：在压力面上，靠近叶片前缘的上半部分的面型对压比比较敏感；在吸力面上，在靠近叶根的部分和叶片中间的面型对压比较敏感。压力面与吸力面相比，吸力面对目标函数更敏感，因此将对目标函数敏感度高的叶片吸力面区域作为优化变形区域，如图3所示。

通过改变叶片敏感区的面型，以达到提高转子总压比的目的。采用一阶最速下降寻优算法，将网格变形中心点的坐标和梯度信息作为寻优算法的输入，最小步长设置为0.000 1，经过20步优化迭代，获得较为理想的结果。文献[14]提到的叶片优化设计方法，要对叶片型线进行参数化，优化方法采用响应面模型，需要18个设计变量，290个样本点，计算耗时，与其相比本文的优化过程只需要几十个迭代步，计算量小，从而缩短了优化设计周期。

图2  目标函数对叶片外形的敏感度分布

Fig. 2  Sensitivity distribution of cost function to blade

图3  叶片变形区域

Fig. 3  Deformation zone of blade

4  优化结果分析

对比了在近最大等熵效率工况点下的叶片优化前后的性能，如表2所示。与原结构相比，优化后的转子叶片的进出口总压比提高0.352%，等熵效率提高0.339%，质量流量提高0.362%。图4所示为在近最大等熵效率工况下，转子叶片出口性能参数沿径向的分布曲线。由图4可知：出口总压比在50%叶高以上有较明显提升，等熵效率在60%叶高以上有提升，在50%叶高以下总压比和等熵效率均没有明显改变。

图4  近最大效率工况下出口气动参数展向分布对比

Fig. 4  Radial plots at outlet near peak efficiency

表2  优化前后性能参数对比

Table 2  Comparison of performance parameters baseline and optimized rotor

图5所示为叶片优化前后的吸力面壁面极限流线。由图5可知：在近最大等熵效率工况下，由于激波和吸力面边界层之间的相互作用，在叶片吸力面中间产生1条分离线，在2条分离线之后出现1条附着线，随后在叶片上部区域靠近尾缘的附近出现附着线，说明叶片流场很复杂。对比优化前后吸力面的极限流线，优化后叶片的极限流线在50%叶高以上出现较为明显的变化，第2条分离线前移，形状也发生变化，靠近尾缘上方的附着线与尾缘的距离变大，形状也发生变化。

图6所示为在近最大等熵效率工况下，不同叶高截面的内部马赫数分布。由图6可见：在30%和50%叶高截面，优化前后叶片通道内流场的激波形态和马赫数分布均没有明显变化，但对于65%，75%和85%叶高截面，优化后的通道内近吸力面的激波前马赫数变大，且激波角度明显变小，原叶片的激波近似垂直于叶片压力面，通道内的激波有从正激波向斜激波转化的趋势，在一定程度上减少了通道内流动损失，提高了转子叶片性能。

图5  近最大效率工况下壁面极限流线对比

Fig. 5  Comparison of limiting streamline near peak efficiency

图6  近最大效率工况下马赫数分布对比

Fig. 6  Comparison of Mach-number distributions near peak efficiency

5  结论

1) 与传统的叶片优化设计方法相比，本文提出叶片优化方法，不需要对叶片型线进行参数化，优化过程计算量小，缩短了优化设计周期。

2) 采用伴随方法，可以获得叶片型面对目标函数的梯度，直观地显示了叶片型面的敏感度区域，有助于准确锁定优化区域。

3) 在近最大等熵效率工况下，优化后转子叶片的等熵效率提高0.339%，总压比提高0.352%，质量流量提高0.362%。优化后叶型能够改变50%叶高以上的激波结构，降低流动损失，对应提高叶片性能。

参考文献：

[1] JAMESON A. Aerodynamic design via control theory[J]. Journal of scientific computing, 1988, 3(3): 233-260.

[2] NADARAJAH S K, JAMESON A, ALONSO J J. Sonic boom reduction using an adjoint method for wing-body configurations in supersonic flow[C]//9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization Conference. Atlanta, GA, 2002: 5547.

[3] KIM S, ALONSO J J, JAMESON A. Multi-element high-lift configuration design optimization using viscous continuous adjoint method[J]. Journal of Aircraft, 2004, 41(5): 1082-1097.

[4] PAPADIMITRIOU D I, GIANNAKOGLOU K C. Total pressure loss minimization in turbomachinery cascades using a new continuous adjoint formulation[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 2007, 221(6): 865-872.

[5] PAPADIMITRIOU D I, GIANNAKOGLOU K C. Compressor blade optimization using a continuous adjoint formulation[C]// ASME Turbo Expo 2006: Power for Land, Sea, and Air. Barcelona, Spain: American Society of Mechanical Engineers, 2006: 1309-1317.

[6] WANG D X, HE L. Adjoint aerodynamic design optimization for blades in multistage turbomachines. Part i: Methodology and verification[J]. Journal of Turbomachinery, 2010, 132(2): 021011.

[7] WANG D X, HE L, LI Y S, et al. Adjoint aerodynamic design optimization for blades in multistage turbomachines. Part ii: validation and application[J]. Journal of Turbomachinery, 2010, 132(2): 021012.

[8] 张朝磊, 厉海涛, 丰镇平. 基于离散伴随方法的透平叶栅气动优化设计[J]. 工程热物理学报, 2012, 33(1): 47-50.

ZHANG Chaolei, LI Haitao, FENG Zhengping. Aerodynamic optimization design of turbomachinery cascade based on discrete adjoint method[J]. Journal of Engineering Thermophsics, 2012, 33(1): 47-50

[9] LI Y, YANG D, FENG Z. Inverse problem in aerodynamic shape design of turbomachinery blades[C]//ASME Turbo Expo 2006: Power for Land, Sea, and Air. Barcelona, Spain: American Society of Mechanical Engineers, 2006: 1535-1543.

[10] RENDALL T C S, ALLEN C B. Unified fluid-structure interpolation and mesh motion using radial basis functions[J]. International Journal for Numerical Methods in Engineering, 2008, 74(10): 1519-1559.

[11] RENDALL T C S, ALLEN C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2009, 228(17): 6231-6249.

[12] RENDALL T C S, ALLEN C B. Reduced surface point selection options for efficient mesh deformation using radial basis functions[J]. Journal of Computational Physics, 2010, 229(8): 2810-2820.

[13] ZILINSKAS A. Practical mathematical optimization: an introduction to basic optimization theory and classical and new gradient-based algorithms by Jan Asnyman[J]. Interfaces, 2005(6): 613-615.

[14] 王祥锋, 韩万金. 跨音压气机转子多目标气动优化设计[J]. 吉林大学学报(工学版), 2010, 40(1): 299-304.

WANG Xiangfeng, HAN Wanjin. Multi-objective aerodynamic optimal design for a transonic compressor rotor[J]. Journal of Jilin University (Engineering and Technology Edition), 2010, 40(1): 299-304.

(编辑  杨幼平)

收稿日期：2015-06-04；修回日期：2015-08-27

基金项目(Foundation item)：国家高技术研究发展计划(863计划)项目(2012AA041304)；高等学校博士学科点专项科研基金资助项目(20110061110022)(Project (2012AA041304) supported by the National High Technology Research and Development Program (863 Program) of China; Project (20110061110022) supported by the Specialized Research Fund for the Doctoral Program of Higher Education)

通信作者：张雷，教授，博士生导师，从事先进加工制造技术研究；E-mail：zhanglei@jlu.edu.cn