刘浩　张雷　李霄琳

(吉林大学 机械科学与工程学院, 吉林 长春 130025)

基于网格变形的叶片多目标气动优化设计*

刘浩张雷†李霄琳

(吉林大学 机械科学与工程学院, 吉林 长春 130025)

摘要：压气机叶片的传统优化设计存在设计变量多和优化周期长等不足.为此,文中利用非均匀有理B样条基函数建立了转子流场网格自由变形参数化方法,并结合改进的拉丁超立方试验设计、Kriging响应面模型及NSGA-Ⅱ多目标遗传算法构建了转子叶片的气动优化设计体系.计算结果表明:在98%的堵塞质量流量工况下,优化后的叶片总压比提高了0.33%,等熵效率提高了0.83%;优化后压气机转子形状为前倾型叶片,提升了转子性能,降低了激波损失;与传统优化设计方法相比,文中优化体系的设计变量明显减少,缩短了优化设计周期.

关键词：压气机;优化设计;网格变形;响应面;多目标遗传算法

随着航空发动机性能要求的不断提升,在现有发动机型号的基础上进行优化改进和缩放是业内普遍采用的方法,优化压气机叶片结构是提升发动机性能的有效途径.计算机技术及计算流体动力学(CFD)技术的发展,使发动机叶片的数值优化成为设计过程中必不可少的环节.国内外学者对此进行了大量研究,Lian等[1]运用遗传算法(GA)、梯度基算法和响应面方法对Rotor67的压比、重量进行多目标优化.Okui等[2]运用人工神经网格遗传算法对压气机叶片进行优化以提高转子效率.李彬等[3]建立了基于自适应差分遗传算法的长叶片透平多学科多目标优化设计系统.

上述叶片优化设计方法是将叶片截面轮廓和积叠线相关的几何变量作参数化,优化迭代的每一步都会重复几何形状生成和网格划分,降低了优化的效率和网格的质量.特别是在复杂几何特征位置处和边界层的网格,容易出现网格划分失败.此外,传统优化方法对几何模型的参数化要求高,几何模型易更新失败.

为解决上述问题,将网格变形[4]引入压气机叶片结构的优化设计中,通过更新网格模型的网格节点位置产生新的形状,不需要对几何模型进行参数化及重复划分网格,提高了计算效率.此外,变形后的网格与原网格具有一致的网格拓扑结构,保证了网格质量,消除了网格差异对计算流体动力学计算结果分析产生的干扰.

目前,运用网格变形技术进行结构优化设计已成为国内外研究的热点.王元元等[5]将自由变形(FFD)技术应用于运输机的气动优化中,数值模拟和风洞试验结果表明,优化结果达到了改善后体流动、减小后体阻力的目的.Yusuke等[6]采用FFD技术对多体船进行多目标优化设计,并对优化后的结果进行了模型试验验证.Biancolini等[7]运用径向基函数网格变形技术对船的帆型进行优化设计,体现了网格变形技术的优越性.

文中利用网格变形技术,结合改进拉丁超立方体试验设计、Kriging响应面模型及多目标遗传算法对NASA Rotor37轴流压气机叶片叶型进行多目标气动优化,以网格变形参数为设计变量,等熵效率、总压比为目标函数进行优化设计.

1数值模型

文中基于Fluent软件平台,以Rotor37为研究对象,对雷诺平均N-S控制方程进行三维定常计算.采用Spalart-Allmaras湍流模型,对流项采用二阶精度的迎风格式进行离散.为保证网格在叶片前、后缘附近有较好的正交性,采用了分区、贴体的网格拓扑结构.流道叶身网格为77×57×209的结构型网格,叶顶间隙采用25×33×209的蝶形网格,网格总数约为155万,并且保证了计算结果的网格无关性.根据试验条件设定边界条件,入口边界给定平均分布的总压、总温;出口边界给定轮毂半径处的静压,其余出口压力按径向平衡方程计算得到;壁面设定为绝热、固壁、无滑移边界条件.

2网格变形及参数化

文中在圆柱坐标系下构造了FFD控制体,控制体的形状与转子流场外形相似,控制体是连续体,控制顶点可以任意分布.

Rotor37流场网格变形的基本步骤如下:

(1)合并3个参数控制体(流场入口至叶片前缘区域、叶片前缘至后缘区域、叶片后缘至出口区域)之间的交界面成为一个控制体.控制体的参数空间为Q(u,v,w)、局部坐标系O-UVW为圆柱坐标系,X0(O)是局部坐标系的原点,U、V、W是轴矢量,则任意一点X在局部坐标系中的坐标表示为

X(u,v,w)=X0+uU+vV+wW

(1)

参数空间坐标为

(2)

(3)

(4)

且满足0

(2)添加控制顶点Pi,j,k,并将待变形网格映射到参数空间,计算每个网格节点X在参数空间中的参数坐标,则任意网格节点在参数空间中的坐标为

(5)

(0≤u,v,w≤1)

(6)

式中, Wi,j,k为Pi,j,k的权因子,Bi,p(u)、Bj,q(v)、Bk,o(w)分别为定义在节点向量上的p阶、q阶、o阶非均匀有理B样条基函数,u={u0,u1,…,up,…,ul,…, ul+p},v={v0,v1,…,vq,…,vj,…,vj+q},w={w0,w1,…,wo,…,wk,…,wk+o}.

(4)保持每个网格节点在参数空间中相应的参数坐标不变,由于控制体的形状发生改变,故每个网格节点的空间位置也会发生改变,从而实现网格的变形.设网格节点X移动到X′,则X′的坐标可表示为

(7)

图1给出了Rotor37流场网格和参数化控制体构造过程.图中的4个控制点接近叶片叶顶部位,改变其位置叶片随之发生明显的变形.相关研究表明,Rotor37的倾斜和掠型叶片有助于提高转子的气动性能[8-10].叶片随4个控制点沿不同方向移动而产生的倾斜和掠型效果如图2所示.表1列出了控制点经过许用最大位移后,得到变形后网格的质量数据,对比发现变形前后的网格质量参数变化很小,变形后的网格与原网格保持同一网格质量精度.

图1　网格变形参数控制体Fig.1　Parametric deformable control volume of mesh

图2　变形后的倾斜和掠型叶片Fig.2　Leaned and swept blades after deformation

表1　原网格和变形后网格的质量对比Table 1　Comparison of quality between the original mesh and the deformed mesh

3优化方法与设计

优化步骤如下:①对所要优化的问题进行数值描述,划分网格;②参数化模型,确立设计变量;③通过试验设计生成样本点;④对样本点进行数值模拟,获得响应值数据库;⑤采用Kriging响应面方法构造近似模型;⑥确立目标函数和约束条件,并采用多目标遗传算法进行全局寻优;⑦对最优设计点进行数值流场验证,判断数值模拟结果与近似模型的响应值偏差是否在误差范围内;⑧将不满足条件的数值模拟解加入数据库,返回步骤⑤直至误差在允许范围内,得到优化结果,最终输出网格文件和几何文件.

3.1　响应面法

选用非线性程度较高和对局部响应突变问题具有较好拟合效果的Kriging模型[11]作为响应面近似模型.该模型已在航空航天、船舶、汽车等工程的优化设计领域得到广泛的应用.

根据Kriging模型原理,可将未知函数关系表示为

y(x)=f(x)+Z(x)

(8)

式中,未知函数f(x)表示全局近似,Z(x)表示服从正态分布N(0,σ2)的随机误差模型,用于局部近似.

随机函数Z(x)的协方差矩阵为

cov[Z(x(i)),Z(x(j))]=σ2R[R(x(i),x(j))]

(9)

式中:i,j∈{1,2,…,ns};R(x(i),x(j))为样本空间中任意两个样本点x(i)和x(j)的空间相关函数;R为ns×ns相关系数矩阵,其元素由相关函数R(x(i),x(j))计算得到,R(x(i),x(j))对模拟精度起关键作用.文中采用高斯函数求解R(x(i),x(j)),即

(10)

(11)

(12)

(13)

3.2　多目标优化算法

文中采用改进的非支配排序方法NSGA-Ⅱ[12]进行多目标优化.多目标优化问题的最优解不再是给定约束条件下使所有目标函数最大或最小,而是Pareto最优解集.Pareto解也叫非支配解,当

(14)

时,解x0支配或优于x1(x0

(15)

Pareto最优前沿指所有Pareto最优解对应的目标函数空间的区域PF,

(16)

NSGA-Ⅱ最突出的特征是采用了非支配排序和适应度共享,前者在优化过程中收敛到Pareto最优前沿,后者则保证了Pareto最优解的多样性.由于文中的目标函数之间存在矛盾,不可能使各个目标函数同时达到最优解,因此优化问题的最终解是从所有Pareto最优解中选择一个最优折衷解.

3.3　设计变量与目标函数

如图3所示,将控制网格变形的控制点作为设计变量,文中按层选取控制点,对应的径向比例系数分别为0.00、0.25、0.60、1.00.将16个控制点作为设计变量,由于网格模型是周期对称的,为保证周期网格交界面在变形后仍能保持一致,必须将控制点成对分组,将靠近叶片前缘的控制点划为组a、c、e、g,叶片尾缘的控制点划为组b、d、f、h.叶片根部控制点组g、h固定不动,控制点组a、b、c、d、e、f分别沿转子叶片的轴向和切向变化以控制网格变形,因此共有12个有效设计变量.与考虑叶型截面及积叠线的传统优化设计方案(文献[13]使用了23个设计变量,文献[14]使用了18个设计变量)相比,文中方案的设计变量明显减少,计算时间相应减少,从而缩短优化设计周期.

图3　控制点的分组Fig.3　Grouping of control points

叶片气动性能的多目标优化设计,以等熵效率η最大为目标函数、总压比r1变化不超过0.5%为约束条件进行优化.

目标函数为

(17)

式中，γ是设计变量,rref是参考总压比.

4结果分析

4.1　数值计算的有效性

利用Denton[15]给出的试验数据校验数值计算的精度和可靠性.湍流模型为Spalart-Allmaras模型,并且保证数值结果的网格无关性,计算堵塞质量流量是20.87kg/s,在试验堵塞质量流量(20.93±0.14)kg/s的偏差范围内.计算和试验得到的转子特性曲线如图4所示，其中r2为流量比.计算得到的总压比在分布趋势和数值上与试验值都有很好的一致性,等熵效率在分布趋势上与试验值一致,但在数值上有差别,峰值等熵效率比试验值低近2%,但这些差别与文献[15]中给出的仿真计算结果基本一致,试验失速点流量比约为0.919 0,数值计算失速点流量比为0.913 5,相差0.59%,因此文中的数值计算结果是可信的.相比其他湍流模型,Spalart-Allmaras(S-A)模型在保证计算精度的同时可获得较高的计算效率,针对优化问题使用Spalart-Allmaras模型是较好的选择.

图4　转子特性的数值计算值和试验值对比Fig.4　Comparison of rotor performance between numerical calculation and experimental values

4.2　优化结果分析

在Rotor37的优化设计中,Benini[13]的优化使转子在特定工况(总压比不变)下的等熵效率提高了1.5%;王祥锋等[14]的优化使转子等熵效率提高了0.7%,总压比提高了0.66%;Samad等[16]的研究使转子等熵效率提高了0.51%,总压比提高了1.25%.

文中采用12个设计变量,试验设计(DOE)方法采用改进的拉丁超立方抽样法,即最优空间填充设计,共生成281个样本点.在98%的堵塞质量流量工况下,采用多目标遗传算法得到的Pareto线如图5所示.图5表明,随着等熵效率的提高,进出口总压比逐渐降低.要确定目标函数最优设计点时,需要先从Pareto分布中挑选折衷的优化候选点,并对每个设计点分别进行CFD计算,再进行分析比较,进而确定最优设计点.在98%的堵塞质量流量计算工况下,等熵效率从优化前的0.863 57变为0.870 73,增加了0.83%,总压比从优化前的2.082变为2.089,增加了0.33%.

图5　多目标遗传算法的Pareto前沿分布Fig.5　Pareto front distribution of multi-objective genetic algorith

图6对比了优化前、后叶片的外形变化,优化后的叶片形状可以看作是前倾型叶片.在不同设计转速(100%、80%、60%)下优化前、后转子的性能曲线如图7所示,从图可以看出,在所有工况下转子的等熵效率都有实质性的提升.以100%设计转速为例,在95%~99%的堵塞质量流量工况范围内,转子效率平均提高近0.7%,而总压比变化不明显.在最大等熵效率工况下出口转子性能参数沿径向的分布曲线如图8所示.优化后转子叶片的等熵效率在20%叶高以上得到提升,其总压比在30%~85%叶高之间有所提高.出口转子性能的变化是由于叶片曲率的改变引起三维激波结构发生了变化,而流场的激波现象对转子的局部性能有主要的影响.

图6　优化前、后叶片形状对比Fig.6　Comparison of blade profile between before and after optimization

图7　优化前、后转子的性能对比Fig.7　Comparison of performance between original and optimized rotor

图8　最大效率工况下出口气动参数展向分布Fig.8　Aerodynamic performance radial plots at outlet at the peak efficiency

图9　最大效率工况下马赫数等值线分布对比Fig.9　Comparison of Mach-number isoline distributions at peak efficiency

在最大等熵效率工况下叶片通道内马赫数等值线分布如图9所示,叶片型面的变化对激波空间形态的改变有明显的作用,而各性能参数及流场结构的变化主要来源于前缘激波的结构变化.在50%叶高处激波并没有明显变化,在70%叶高处激波变为后掠激波,相应的激波损失减小,流道内流体马赫数从原来的0.9增加到1.1;在90%叶高处激波是后掠激波,前缘激波无论在强度还是在波形上都发生了明显变化,流道内流体马赫数由1.0增加到1.2,由强正激波变为弱斜激波,阻力减少,因此提升了转子效率.

图10显示了叶片优化前、后的壁面极限流线.原叶片因激波和吸力面边界层之间的相互作用而在叶片型面接近中间处产生一条分离线,在边界层分离区域靠近叶根的部位分离线之后出现一条附着线,在附着线之后又出现了分离线,随后在叶片上部区域靠近尾缘的附近出现附着线,说明原叶片下游流场复杂,出现二次分离,边界层分离严重,分离泡多,损失大.优化后叶片的壁面极限流线与原叶片相比,分离线在50%叶高以上开始明显向吸力面下游倾斜,说明激波由正激波逐渐变为斜激波;分离区域变小,并且在叶片上部靠近尾缘的区域没有出现附着线,说明分离区域内的流场损失比原叶片小,因此优化后叶片性能得到了提高.

图10　最大效率工况下壁面极限流线对比Fig.10　Comparison of the wall surface limiting streamline at the peak efficiency

5结论

文中针对NASARotor37叶片的优化获得以下结果:

(1)运用网格变形技术结合响应面模型及多目标优化算法对压气机转子叶片进行气动优化,叶片性能均得到一定的提升;

(2)基于自由变形控制体的网格变形方法的计算效率高,变形后的网格质量好,节省了重新划分网格的时间,同时相比叶型截面及积叠线的优化设计方法,输入变量的数量明显减少,缩短了优化设计周期;

(3)在98%的堵塞质量流量工况下,优化后转子的等熵效率提高了0.83%,总压比提高了0.33%,优化后的叶型为前倾叶片,结果表明前倾型叶片能显著改善叶片前缘的激波结构,削弱激波强度,降低激波损失,提升转子性能.

参考文献:

[1]LianY,LiouMS.Multi-objectiveoptimizationoftranso-niccompressorbladeusingevolutionaryalgorithm[J].JournalofPropulsionandPower,2005,21(6):979-987.

[2]OkuiH,VerstraeteT,vandenBraembusscheRA,etal.Three-dimensionaldesignandoptimizationofatransonicrotorinaxialflowcompressors[J].JournalofTurbomachinery,2013,135(3):031009.

[3]李彬,宋立明,李军,等.长叶片透平级多学科多目标优化设计 [J].西安交通大学学报,2014,48(1):1- 6.

LiBin,SongLi-ming,LiJun,etal.Multidisciplinaryandmulti-objectiveoptimizationdesignoflongbladeturbinestage[J].JournalofXi’anJiaotongUniversity,2014,48(1):1- 6.

[4]张湘玉,廖文和,刘浩.保细节的基于曲面控制的网格变形术 [J].华南理工大学学报:自然科学版,2010, 38(7):39- 44.

ZhangXiang-yu,LiaoWen-he,LiuHao.Detailpreservingmeshdeformationbasedonsurfacecontrol[J].JournalofSouthChinaUniversityofTechnology:NaturalScienceEdition,2010,38(7):39- 44.

[5]王元元,张彬乾,郭兆电.基于FFD技术的大型运输机上翘后体气动优化设计 [J].航空学报,2013,34(8):1806-1814.

WangYuan-yuan,ZhangBin-qian,GuoZhao-dian.AerodynamicoptimizationdesignforlargeupsweptafterbodyoftransportaircraftbasedonFFDtechnology[J].ActaAeronauticaetAstronautiaSinica,2013,34(8):1806-1814.

[6]YusukeT,DanieleP,EmilioFC,etal.Singleandmultiobjectivedesignoptimizationofafastmultihullship:numericalandexperimentalresults[J].JournalofMarineScienceandTechnology,2011,16(4):412- 433.

[7]BiancoliniME,ViolaIM,RiotteM.SailstrimoptimisationusingCFDandRBFmeshmorphing[J].Computers&Fluids,2014,93(10):46-60.

[8]BeniniE,BiolloR.Aerodynamicsofsweptandleanedtransoniccompressor-rotors[J].AppliedEnergy,2007, 84(10):1012-1027.

[9]JangCM,SamadA,KimKY.Optimaldesignofswept,leanedandskewedbladesinatransonicaxialcompressor[C]∥ASMETurboExpo2006:PowerforLand,Sea,andAir.Barcelona：AmericanSocietyofMechanicalEngineers,2006:1279-1288.

[10]BiolloR,BeniniE.Recentadvancesintransonicaxialcompressoraerodynamics[J].ProgressinAerospaceSciences,2013,56:1-18.

[11]SimpsonTW,MaueryTM,KorteJJ,etal.Krigingmodelsforglobalapproximationinsimulation-basedmultidisciplinarydesignoptimization[J].AIAAJournal, 2001,39(12):2233-2241.

[12]DebK,PratapA,AgarwalS,etal.Afastandelitistmultiobjectivegeneticalgorithm:NSGA-II[J].IEEETransactionsonEvolutionaryComputation,2002,6(2):182-197.

[13]BeniniE.Three-dimensionalmulti-objectivedesignoptimizationofatransoniccompressorrotor[J].JournalofPropulsionandPower,2004,20(3):559-565.

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

WangXiang-feng,HanWan-jin.Multi-objectiveaerodynamicoptimaldesignforatransoniccompressorrotor[J].JournalofJilinUniversity:EngineeringandTechnologyEdition,2010,40(1):299-304.

[15]DentonJD.Lessonsfromrotor37 [J].JournalofThermalScience,1997,6(1):1-13.

[16]SamadA,KimKY.Shapeoptimizationofanaxialcompressorbladebymulti-objectivegeneticalgorithm[J].ProceedingsoftheInstitutionofMechanicalEngineers,PartA:JournalofPowerandEnergy,2008,222(6):599- 611.

Multi-Objective Aerodynamic Optimum Design of a Blade Based on Mesh Deformation

LiuHaoZhangLeiLiXiao-lin

(School of Mechanical Science and Engineering, Jilin University, Changchun 130025, Jilin, China)

Abstract:There exist the defects of multi design variables and long optimization cycle in the traditional optimum design of compressor blades. In order to solve these problems, this paper proposes a free-form mesh deformation parameterization method of fluid grids by using a non-uniform rational B-spline basis function, and constructs an design system of aerodynamic optimization of rotor blades by combining the advanced design of Latin hypercube sampling experiments, the Kriging response surface model and the NSGA-Ⅱ multi-objective genetic algorithm. Calculation results show that (1) the optimized blade has a total pressure ratio improvement by 0.33% and an isentropic efficiency improvement by 0.83% at a choke mass flow of 98%; (2) the optimized blade is a forward-leaned blade, which helps reduce the shock loss and improve the performance of the rotor; and (3) in comparison with the traditional optimization design method, the proposed optimization system reduces design variables and shortens optimized cycles.

Key words:compressor; optimum design; mesh deformation; response surface; multi-objective genetic algorithm

中图分类号：V232.4

doi：10.3969/j.issn.1000-565X.2015.09.022

作者简介：刘浩(1985-),男,博士生,主要从事叶片气动优化设计研究.E-mail: adj1212@163.com† 通信作者： 张雷(1970-),男,教授,博士生导师,主要从事先进加工制造技术研究.E-mail: zhanglei@jlu.edu.cn

*基金项目：国家“863计划”项目(2012AA041304)；高等学校博士学科点专项科研基金资助项(20110061110022)

收稿日期：2014-10-28

文章编号：1000-565X(2015)09-0141-08

Foundation items： Supported by the National High-Tech R & D Program of China(2012AA041304) and the Specialized Research Fund for the Doctoral Program of Higher Education(20110061110022)