基于伴随方法的单级低速压气机气动设计优化

2020-06-03罗佳奇杨婧

航空学报 2020年5期

罗佳奇，杨婧

1. 浙江大学航空航天学院，杭州 310027 2. 德州农工大学机械工程系, College Station TX 77843

随着高性能计算的快速发展，基于计算流体力学(CFD)的优化设计技术在叶轮机械叶片气动设计中已经得到了应用。基于CFD的气动设计优化(ADO)不仅能够缩短研制周期、降低成本，相对于依靠设计者经验的设计方法还能提高设计结果精度。对于压气机而言，通过优化叶片气动外形增加压比、降低流动损失对于提高压气机气动性能具有重要意义。

目前得到应用的优化方法分为两类：梯度方法和非梯度方法。进化算法[1-3]和代理模型[4]是典型的非梯度方法。Benini[2]采用进化算法对跨声速压气机转子叶片NASA Rotor 37的中弧线和厚度分布进行了优化设计，Lian和Liou[3]采用进化算法和代理模型对跨声速压气机转子叶片NASA Rotor 67进行了多目标优化设计，Samad等[4]采用多层次代理模型对转子叶片的叶片掠、叶片倾斜及叶片弯进行了优化设计。上述非梯度方法虽然在叶轮机械ADO中已经得到应用，但是对于多设计参数的复杂ADO问题，尤其是需要通过求解Navier-Stokes方程计算流场时，将需要大量的计算资源，制约了这类方法的工程应用。

基于灵敏度分析的优化方法在叶轮机械ADO中得到了应用，直接差分法(FDM)和伴随方法是两种常用的灵敏度计算方法。采用FDM计算灵敏度时，所需的流场计算次数和设计参数线性相关。在20世纪80年代，Jameson[5]成功地将伴随方法应用于飞行器气动外形优化设计。采用伴随方法，每个气动函数的灵敏度计算只需要大约两次流场计算，计算时间和设计参数个数基本无关。之后，基于伴随方法的翼型、机翼以及翼身组合体ADO研究得到了迅速发展[5-9]。由于灵敏度计算的高效性和精确性，伴随方法在叶轮机械ADO中的应用也越来越广泛[10-19]。Luo等[17-18]开展了基于伴随方法的跨声速压气机叶片的多工况、多目标ADO研究。近年来，掺混面模型方法也应用于伴随方程的数值求解，推动了基于伴随方法的叶轮机械多排ADO研究的发展[13,14,16,20]。

本文将发展一种基于掺混面模型的叶轮机械多排伴随方法，并应用于某型低速、低压比压气机ADO。首先采用自主开发的初步设计方法设计4.5级压气机叶片的气动外形，通过数值模拟发现：最后级静子叶片的端壁和吸力面角区内存在较大分离，严重影响了压气机效率。之后介绍适用于多排ADO的伴随方法，并在近失速工况下对最后级静子叶片的中弧线及安装角进行优化，优化目标是降低出口熵增，同时约束出口流量。最后开展多工况ADO，在两个不同转速下进行最后级静子叶片的外形优化设计，分析外形优化设计对设计工况、非设计工况气动性能的影响。

1 压气机初步设计

该低速4.5级压气机由第1排进口导叶和后4级叶排构成，给定设计压比和等熵效率。在初步设计之前需要确定一些前提条件：单级多变效率和压气机相等，给定单级温比，压气机进口到导叶出口的面积收缩率等。采用基于经验修正[21]的4.5级压气机初步设计的基本流程为：

第1步给定效率初值，确定压气机进、出口流道几何参数、速度三角形和气动参数。

第2步根据给定的温升分布、静子出口气流角和流道参数等，进行逐级设计，迭代计算各级转子、静子出口截面的相关参数。

第3步计算整机压比和效率，判断是否满足设计条件，如果不满足则用计算效率替代初始效率，重复设计。

第4步采用简单径向平衡方程确定轮毂和机匣的相关参数。

第5步确定压气机轴向尺寸，进行二维叶片造型和径向堆积。

图1为该4.5级压气机轮毂、叶中和机匣处的叶型，图中：x表示轴向，Rθ表示周向。值得注意的是：初步设计时，导叶和静子叶片均采用直叶片设计。流场分析发现：直叶片设计使得在前三级的影响下，最后级静子叶片的端壁和吸力面角区流动分离严重；因而，该初步设计的4.5级压气机存在较大优化空间，可以通过叶片外形优化设计抑制流动分离，提高压气机气动性能。

图1 3个不同径向位置4.5级压气机叶片叶型Fig.1 Blade profiles of the 4.5-stage axial compressor at three different spans

2 气动设计优化

ADO的流场计算采用流场求解器Turbo90[22-24]，给定目标函数及优化约束的数学建模以确定伴随方程的边界条件并实现数值求解。确定伴随灵敏度后，采用简单、易实现的最速下降法来进行优化设计，气动外形扰动步长根据经验给定。

2.1 流场数值模拟

压气机最后级的流场通过求解雷诺平均Navier-Stokes (RANS)方程及Spalart-Allmaras (SA) 湍流模型[25]确定。采用多重网格、隐式残值光顺、当地时间步长等技术加速计算收敛。最后级进口流动参数如总压、总温、流动角由4.5级压气机流场数值模拟结果给定，出口处给定轮毂压力并由径向平衡方程及其他流动参数确定压力径向分布，转-静交接面上采用掺混面模型[26]实现流场连续求解，并采用无反射边界条件[27]。本文所采用的Mixed-out掺混面模型能够保证交接面上流量、动量和能量的通量守恒，其详细数学形式参考附录A。

计算采用H-型结构网格，不考虑叶尖间隙影响，S1及S2流面网格如图2所示。首先在3套不同的网格上进行流向和周向网格无关性研究，径向网格数固定为73，由此确定收敛的流向网格数；之后在3套不同网格上进行周向和径向网格无关性研究，叶面流向网格数固定，周向和径向网格增长比为2。网格无关性研究的具体实现流程可参考本文作者前期工作[19]。表1为近失速工况下周向和径向网格无关性研究结果，包括：归一化流量、静子效率、级效率、压比。由表可知：相对于网格2与网格1的流场偏差，网格3与网格2的流场偏差明显减小，流量、级效率的相对偏差分别为0.31%和0.46%，压比偏差接近于零。后续研究将采用网格2，转子和静子的网格规模均为：叶面流向97、周向49、径向57个网格点。

实际上，为了进一步证实网格无关性解的有效性，研究中对优化后的最后级也采用了相同规模的密网格(网格3)进行流场数值模拟，其气动性能也得到了提升。

图2 最后级网格Fig.2 Grids for the last stage

表1 网格无关性解Table 1 Independent grid solutions

2.2 参数化建模

中弧线、厚度分布、安装角是叶片的3个重要设计参数。本文的ADO研究中，只改变最后级静子叶片的中弧线和安装角径向分布，其厚度分布不变。中弧线的变化通过叠加Hicks-Henne型函数[28]实现，安装角变化在前期工作中已经介绍[15]。选取10%、15%、20%、30%、70%、80%、85%、90%这8个径向位置，每个径向位置的安装角都作为设计参数，叶型中弧线变化由8个型函数叠加确定；共有72个设计参数。

图3为叶型中弧线变化和安装角变化的几何示意图，其中：δ表示变化量，下标 camber 和 stagger 分别表示中弧线和安装角。

叶型变化后，在原始网格基础上采用插值网格扰动技术[7]重新生成计算网格：

xnew=xold+r(xW,new-xW,old)

(1)

式中：x表示网格坐标；下标 old和 new 分别表示原始和更新后的网格,W表示叶面网格点；r为归一化的空间距离，在叶面网格点r=1。

图3 设计参数变化Fig.3 Design parameter variations

2.3 目标函数

本文ADO的主要目标是降低流动损失，同时约束出口流量。叶轮机械的流动损失可以用熵增来衡量，因而目标函数定义为

I=Sgen+Λ|σ-1|

(2)

式中：σ表示设计流量与原始流量的比值；Λ为约束系数；Sgen表示平均熵增，定义为

(3)

其中：Δs表示每个单元上的熵增，其定义见文献[15]。

基于伴随方法的多工况ADO方法在前期工作中已经介绍[17]，定义加权目标函数为

(4)

式中：I表示气动参数；M为工况数；λi表示权重且满足 ∑λi=1。相应地，灵敏度计算表达式为

(5)

式中：Gi表示不同工况下的气动参数灵敏度。

本文研究采用最速下降法作为优化方法，所需要的灵敏度由伴随方法计算。接下来将重点介绍伴随方法的基本原理。

3 伴随方法

伴随方法基本原理较为成熟，这里只作简单介绍。定义I为气动参数，一般与流场w和气动外形相关，而气动外形与设计参数V相关，即

I=I(w,V)

(6)

气动参数关于设计参数的灵敏度为

(7)

式中：右端第1、第2项分别表示流场变分、气动外形变分对灵敏度的影响，灵敏度计算的关键是流场变分的计算。如果采用FDM计算灵敏度，虽然原理简单，但是计算效率偏低。值得注意的是，流场和设计参数满足流体力学方程R，即

R(w,V)=0

(8)

类似于方程(7)，方程(8)可以写成

(9)

伴随方法的主要目标是通过将流体力学方程作为约束引入到灵敏度计算中，消除流场变分对灵敏度计算的影响。定义伴随变量Ψ=[Ψ1,Ψ2,Ψ3,Ψ4,Ψ5]T，将伴随变量与方程(9)相乘后再与方程(7)相减，有

(10)

若方程(10)右端第1项的系数为零，即

(11)

则有

(12)

方程(11)称为方程(8)的伴随方程，Ψ为伴随变量，方程(12)为灵敏度计算公式且与流场变分无关。由此可知：只需求解一次流体力学方程和一次伴随方程，即可由方程(12)计算气动参数I的灵敏度。此外，伴随方程是线性方程组，可以采用和RANS方程相同的数值格式进行求解[11,29]。

Giles和Pierce[29]指出了伴随方程进、出口边界条件确定方法：伴随方程的特征线方向与流体力学方程相反。对于多排压气机叶片，转子出口、静子入口的伴随方程边界条件需要采用掺混面方法[13,16]确定。本文采用基于Mixed-out掺混面模型的伴随掺混面方法，首先计算伴随通量Fa1、Fa2、Fa3、Fa4和Fa5：

(13)

式中：r、p、H分别表示密度、压力和总焓；ds为面积单元；u=[ux,uθ,ur] 表示速度矢量；n为单位法向矢量。之后由方程(13)可以计算平均伴随变量Ψa，即

(14)

转-静交接面上仍采用Giles和Pierce提出的无反射边界条件。具体实现流程为[10,13]：在亚声速交接面下游(对应下游静子入口)，伴随特征变量Ψa1、Ψa2、Ψa3、Ψa4由内场确定，Ψa5由交接面上游确定；在亚声速交接面上游(对应上游转子出口)，Ψa1、Ψa2、Ψa3、Ψa4由交接面下游确定，Ψa5由内场外插确定。确定所有伴随特征变量后，伴随变量Ψa1、Ψa2、Ψa3、Ψa4、Ψa5可通过求解线性方程组确定。

4 结果分析

在基于伴随方法的ADO研究中，首先需要对伴随灵敏度进行精度验证，将基于FDM的灵敏度收敛解作为精确解。图4给出了30%叶高的灵敏度分布，前8个为中弧线设计参数灵敏度，第9个为安装角变化灵敏度；AD表示伴随方法。

整体上，从叶片前缘(LE)至叶片中部，伴随灵敏度与FDM灵敏度较为接近；从叶片中部至叶片尾缘(TE)，伴随灵敏度与FDM存在一定偏差，主要原因是没有考虑湍流模型方程的伴随方程求解，而是采用“定涡”假设(Frozen Eddy Viscosity)，即忽略湍流变量变分对灵敏度的影响。实际上，当流动分离较大时，忽略湍流影响会降低伴随灵敏度精度[30]。图4中，后3个中弧线设计参数(6,7,8)正处于静子叶片分离区。虽然精度有所下降，但是伴随灵敏度的分布与FDM基本一致，本文将采用基于“定涡”假设的伴随方法进行ADO研究。

图4 灵敏度对比Fig.4 Comparisons of gradients

4.1 近失速点单点优化

该4.5级压气机最后级静子角区存在流动分离，且在接近失速工况时分离区愈大。单点优化设计在设计转速近失速工况下进行，约束函数系数Λ= 0.4。图5给出了等熵效率和流量的收敛曲线，流量最大相对偏差不超过0.1%，约束效果明显。表2给出了优化前后的熵增、流量及静子和最后级的等熵效率，优化后静子和级效率分别增加了约0.29%、0.15%。

研究中采用相同规模的密网格(表1中的网格3)对优化叶片流场进行数值模拟，表3给出了优化前后密网格的计算结果。与原始叶片相比，优化后静子和级效率分别增加了0.28%、0.13%。这种较好的网格一致性结果进一步证实了本文网格无关性验证的有效性及优化设计结果的可靠性。

图6为优化前后静子叶片的气动外形。由图6(a) 可知：优化设计后自前缘至20%弦长，叶片中弧线曲率增加；自20%弦长至90%弦长，叶片中弧线曲率降低；叶片气动外形的这种变化将使前缘附近流动加速、来流角与叶片前缘角偏差降低，有利于抑制叶背流动分离。

图5 单点优化熵增和流量收敛曲线Fig.5 Convergence history of entropy and mass flow rate under one point optimization

表2 不同叶片的气动性能参数Table 2 Aerodynamic performance parameters of different blades

表3 密网格气动性能参数Table 3 Aerodynamic performance parameters of fine grid

为了比较分析优化设计对最后级流场的影响，图7给出了静子吸力面附近的流线，图8给出了最后级10%叶高S1流面的马赫数云图。由图7可知：原始静子轮毂、机匣角区均存在流动分离，轮毂角区尤为明显；优化设计使轮毂流动分离区朝尾缘和端壁移动，流动分离得到了较好抑制；机匣流动分离区也得到了较好抑制。由图8可以更好地看出优化设计对静子叶片叶背流动分离的抑制作用。整体上，对于低速压气机，流动分离区的减小将是效率增加的主要原因。

图6 优化前后静子叶片气动外形Fig.6 Stator blade aerodynamic shape before and after optimization

图7 静子吸力面附近流线Fig.7 Streamlines near stator suction side

上述优化设计能够改善最后级在指定工况下的气动性能，在压气机全工况范围内对气动性能的影响更值得关注。图9给出了最后级、静子叶片的等熵效率η特征工作线。虽然在其他工况条件下气动性能改善并不明显，但是静子、级效率在全工况范围内都得到一定的提升。若对最后级的转子叶片也进行优化设计，甚至对整个4.5级压气机进行优化设计，通过调整上游叶排对下游流动的影响，将能进一步提升压气机气动性能。

值得注意的是：单点优化设计的工况点远离失速边界，叶背流动分离区较小，通过控制流动分离获得效率提升有限。实际上，在前期工作中采用基于本征正交分解的混合模型方法对该型号压气机最后级也进行了优化设计[31]：工况点距离失速边界较近，流动分离较剧烈，通过最后级静子安装角优化使等熵效率提高了1.47%。因此，虽然通过本文的优化设计获得的效率提升较为有限，但是由上述分析可知：优化设计结果是可靠的，在工程中仍具有应用价值。

图8 10%叶高S1流面马赫数云图Fig.8 Mach number contours at 10% span S1 stream surface

图9 优化前后效率工作线Fig.9 Adiabatic efficiency for both original and optimized blades

4.2 多工况优化

多工况优化设计研究中，将对最后级静子叶片在100%设计转速的近失速工况(P1)、140%设计转速的近失速工况(P2)进行中弧线和安装角的优化设计。在每个转速下分别计算其灵敏度，最后由方程(5)计算加权灵敏度，并由最速下降法实现优化设计,权重λ1=λ2=0.5。

图10给出了两个不同转速下(P1和P2)熵增和出口流量的收敛曲线。表4给出了优化前后叶片气动性能。值得注意的是：表4中 P2 转速下原始和优化叶片的熵增和流量都统一采用 P1 转速下原始叶片的熵增和流量进行归一化处理。多工况优化设计使静子效率、级效率在 P1 转速下分别提升了0.29%和0.15%，在 P2 转速下分别提升了0.24%和0.13%，而两个转速下流量偏差都低于0.15%，表明两个转速条件下流量都得了较好约束。此外，由表4可知高转速条件下熵增增加，其主要原因是：高转速条件下压气机负荷增大，角区分离加剧，流动损失增大。

图10 多工况优化熵增和流量收敛曲线Fig.10 Convergence history of entropy and mass flow rate under multi-points optimization

图11给出了P2转速下最后级10%叶高S1流面的马赫数云图。值得注意的是：图11的最大相对马赫数为0.45，远高于图8的0.375。对比图11和图8发现：随着转速的提高，原始叶片叶背分离更严重，如图8(a)和图11(a)所示；高转速下，通过优化设计静子吸力面的分离区也得到了较好抑制；此外，相对于图8，高转速下分离区的控制效果更加明显。

上述结果表明：通过单点、多工况气动外形优化设计能够提升近失速工况最后级的等熵效率。非设计转速下气动外形设计优化对设计转速下气动性能的影响需要进一步深入分析。图12给出了单点、多工况优化后静子、最后级在非设计转速全工况范围内的效率变化。整体上，多工况优化设计能使静子、最后级的等熵效率在非设计转速全工况范围内得到提升。主要原因是：静子端壁和吸力面角区的流动分离在不同转速条件下都存在，通过气动外形设计优化抑制流动分离的基本原理是一致的：减小来流角和叶片前缘角的偏差、流动加速等。此外，由于ADO在近失速工况下进行，当工作点远离失速边界时，效率提升不断降低；当工作点向阻塞边界靠近时，效率提升不断增加，这一点有待深入研究。