考虑短舱安装参数的全机气动外形优化设计

2019-04-04许丹阳白俊强雷锐午

西安交通大学学报 2019年4期

许丹阳,白俊强,雷锐午

(西北工业大学航空学院,710072,西安)

近年来,高性能计算机硬件平台发展突飞猛进,基于雷诺平均Navier-Stokes方程的气动优化设计方法在飞行器外形设计领域得到越来越广泛的应用[1-3]。

宽体客机作为当前民航市场上的最新一代客机,巡航马赫数达到0.85,处于典型的跨声速非线性区间。客机机翼外形的微小改变往往会产生较大的气动力变化;同时全机多部件之间的互相干扰,也使得气动外形设计工程师无法进行精确定量的外形设计。为了能够合理地对设计空间进行探索,以得到气动特性较好的构型,跨声速客机气动外形优化设计问题往往需要大规模设计变量[4];同时也需要考虑多部件之间的耦合来获得综合性能最优的气动外形。

国外大量的研究成果表明,基于伴随的气动优化设计方法在复杂大规模气动外形优化设计问题中有较好的应用前景。Alonso等基于斯坦福大学非结构网格求解器(SU2)实现了旋翼以及全机构型的气动优化设计[5];Zingg等采用离散伴随方程结合Newton-Krylov隐式求解方法,开展了气动外形优化设计研究取得了显著的成果,并将该方法拓展到了自然层流翼型设计中,实现了层流翼型的优化减阻[6];Martins等利用高可信度多学科优化设计平台(MACH)进行了机翼、翼身组合体以及风力机的优化设计研究[7]。

国内也有学者进行了大量的研究工作,其中左英桃等利用离散共轭方法对DLR-F6机翼、机身、短舱、挂架构型的机翼和短舱进行了几何外形参数化与气动外形优化设计[8];杨旭东等进行了基于黏性伴随方法和Navier-Stokes方程的跨声速机翼气动优化设计研究[9];陈颂等应用伴随方法结合RANS方程进行了民用客机机翼、机身、平尾构型的俯仰力矩配平气动外形优化设计[10];卢娟和朱玉杰分别利用伴随方法对叶栅进行了熵增优化和表面压力分布反设计优化,均有效减小了激波和分离强度,验证了伴随方法在复杂流动中的应用前景[11-12]。

本文在国内学者研究的基础上,选择大型客机短舱下沉量、前伸量、内偏角作为短舱安装参数变量,针对全机构型进行考虑多部件耦合以及俯仰力矩配平的气动外形优化设计,研究了短舱安装参数变化对机翼气动特性的影响。

1 流场数值模拟

在进行三维气动外形优化设计过程中,流场求解采用了Navier-Stokes方程作为流场控制方程。

在笛卡尔坐标系(x1,x2,x3)下,定义速度分量为(u1,u2,u3),采用求和约定惯例,不计热源的三维Navier-Stokes方程守恒形式为

(1)

式中:W是守恒变矢量,Fi为无黏矢通量,Fvi为黏性矢通量,具体表达式如下

(2)

求解Navier-Stokes方程采用自行研制的求解程序[13],其空间离散采用二阶中心差分格式的有限体积法,采用LU-SGS方法进行隐式求解,湍流模型为S-A(Spalart-Allmaras)模型。在求解过程中采用了当地时间步长、多重网格、残值光顺等手段加速流场收敛,程序求解精度参考文献[13]。

2 离散伴随方程法

对于气动外形优化设计问题,一般可以通过对计算网格上的流场变量进行数值积分,将目标函数,例如升力、阻力、力矩表示为如下形式

I=I(W(x),G(x))

(3)

式中:x为每一个优化设计问题所采用的设计变量,取决于所选取的参数化方法;G为所有网格点坐标构成的向量。对于一组给定的设计变量x,通过参数化方法可以得到气动外形对应的表面网格,然后在表面网格的基础上结合动网格方法可以得到空间网格G(x),在其基础上进行流场求解可以得到流场解向量W(x),最后对每个网格上的流场解积分可以得到目标函数I。

当对流场控制方程式(1)进行定常流场求解时,其收敛的流场解残差可以写为

R(W(x),G(x))=0

(4)

将式(3)(4)分别对x求全导数,可以得到

(5)

(6)

为了获得目标函数对于设计变量的导数,需要对式(5)进行求解。当设计变量较多时,为了获得式(5)中dW/dx项,需要针对每个设计变量x进行流场扰动求解。为了节省计算量,可以将式(6)作为约束条件引入式(5),消除dW/dx实现流场计算次数与设计变量个数解耦,得到

(7)

设Φ为伴随算子,整理得到

(8)

最终式(7)可写为

(9)

通过式(8)引入伴随算子,将大规模矩阵求逆问题转换为大规模线性方程组求解问题。式(7)中∂I/∂G为目标函数对空间网格的导数,在收敛流场解的基础上进行网格扰动求解;∂R/∂G为流场残差对空间网格的导数,由数值模拟求解器输出得到;dG/dx为空间网格对设计变量的导数,通过链式法则分解为dG/dS和dS/dx,其中S为参数化方法所确定的气动外形表面网格矢量,dG/dS由动网格方法插值矩阵得到,dS/dx由参数化方法得到。因此,能否高效精确获取目标函数对于设计变量的导数,关键在于大规模线性方程组式(8)的求解。本文采用开源工具箱PETSC中的广义最小残量法(GMRES)[14]来求解伴随方程组(8)。

3 FFD参数化方法

考虑短舱安装参数以及平尾配平影响的优化设计需要同时对多个部件进行参数化,然后通过几何参数化方法定义几何设计变量实现对几何外形的扰动。多块拼接自由变形(FFD)参数化方法通过参数化几何变形量而非几何外形本身,实现对任意外形的参数化,同时由于其光滑、连续、易于求导等特点,在基于梯度的气动外形优化设计问题中应用较为广泛[15],因此本文选取FFD参数化方法。

FFD参数化方法通过下式建立参数空间与物理空间的映射关系[15]

(10)

式中:Bil(r)、Bjm(s)和Bkn(t)分别为l、m和n次Bernstein多项式基函数;Ci,j,k为控制框的控制点在物理空间内的坐标;r、s、t为控制框内的局部坐标;D(x,y,z)为初始表面网格点物理坐标。

Bil(s)的表达式如下

(11)

若已知几何外形在控制框内的局部坐标和变形后控制点在物理空间的坐标,就可以通过下式获得几何变形后在物理空间内的坐标

D′(x,y,z)=

(12)

式中:D′(x,y,z)为变形之后的表面网格点坐标;ΔCi,j,k为FFD控制点的变形量。

针对本文研究的全机构型,参数化控制框如图1所示,全机多部件参数化共采用了43个FFD控制体,相邻的控制体采用点对点对接形式进行拼接,以保证控制体交界面处的几何变形是连续的。43个控制体中3个控制体完成机翼的参数化,一共有462个控制点,可以对机翼外形进行扰动,其设计变量为这些FFD控制点的Z向位移,坐标系如图1所示。机翼FFD控制体中控制点分布在沿展向的11个控制面上,每个控制面上下沿弦线方向各有21个控制点,分别对应上下表面的外形扰动。11个控制剖面分别绕前缘线扭转实现机翼控制剖面扭转角的变化。由于跨声速机翼的压力分布形态、激波位置等对于机翼前缘外形非常敏感,机翼前缘的微小变化会导致气动性能的显著变化,所以FFD控制点在前缘附近进行了适当的加密处理。

图1 全机FFD控制点分布

(a)平尾正偏5°表面网格

平尾FFD控制体绕平尾转轴整体旋转实现平尾的偏转,以配平全机的俯仰力矩。为了保证优化过程中平尾与机身尾部较好地衔接,在平尾与机身相交处添加一个控制体,并且该控制体随平尾的转动而转动,从而保证平尾与机身相交面在优化过程中实现光滑过渡。图2给出了通过平尾控制体旋转得到的平尾偏转构型表面网格,表明在优化过程中通过FFD控制框的整体偏转来实现平尾力矩配平的方案是可行的。

(b)平尾负偏5°表面网格图2 平尾偏转正(上)负(下)5°表面网格

图3给出了短舱局部坐标轴和安装参数示意图。短舱参数化控制框(图1中的机翼下方控制体)整体沿着全机坐标系X方向的前后移动或沿Z方向的上下移动可以对短舱前伸量或短舱下沉量进行扰动,如图3b所示。控制框绕着短舱坐标系的转动可以对短舱安装角进行扰动。短舱坐标系如图3a所示,绕短舱坐标轴Y轴转动可以对短舱内偏角进行扰动,如图3c所示。图4给出了通过偏转短舱FFD控制体实现短舱内偏角变化构型与初始构型的表面网格的对比,从图4可以看出,短舱表面网格随控制框的变化发生了连续偏转,说明本文的参数化方法能够对短舱的内偏角实现参数化扰动。

(a)短舱局部坐标系

(b)短舱前伸量和下沉量

(c)短舱内偏角图3 短舱局部坐标轴和安装参数示意图

图4 短舱内偏角变化构型与初始构型表面网格的对比

4 优化设计流程

本文所采用的优化设计系统主要由多块拼接自由变形(FFD)参数化模块、距离倒数插值(IDW)动网格模块[16]、RANS方程求解计算模块、伴随方程求解模块以及序列二次规划(SQP)优化算法模块[17]构建,优化设计步骤如下,流程如图5所示。

图5 优化设计流程图

步骤1优化算法产生一组设计变量(优化开始时产生的设计变量为0)传递给几何参数化模块,通过几何参数化模块中的FFD参数化方法对需要优化的几何外形进行扰动,产生新的气动外形表面网格。

步骤2动网格模块在新的表面网格基础上利用IDW网格变形算法对空间网格进行插值,得到新的空间网格用于气动求解器的CFD评估分析。

步骤3气动求解器对新的空间网格进行流场分析之后,得到每个网格点上的流场解以及所需要的气动力和力矩特性。

步骤4伴随方程求解模块在收敛的流场解基础上进行伴随方程的构造,利用PETSC中的GMRES算法进行伴随方程的求解,然后通过式(9)获得目标函数对于设计变量的导数。

步骤5将气动目标函数和设计变量导数传递给优化算法,如果优化收敛,则停止优化,输出优化结果,否则跳转至步骤1重复执行。

5 全机气动外形优化设计

5.1 初始构型几何

采用的全机构型(baseline)如图6所示,俯抑力矩参考点位于机翼平均气动弦长(MAC)25%的位置,反映了飞机重心相对于全机的位置,坐标为(30.0,0.0,0.0) m。全机半模参考面积为197 m2,平均气动弦长为8.2 m,机翼前缘后掠角为34.5°,半展长为30 m。宽体客机的设计点巡航马赫数为0.85,升力系数为0.48,基于平均气动弦长的飞行雷诺数为4.9×107。

图6 全机构型初始外形

5.2 不考虑短舱位置变化的优化设计

首先针对宽体客机全机构型进行考虑俯仰力矩配平的减阻优化设计,并将该优化算例命名为优化算例1,其优化模型见表1。

优化算例1经过30轮迭代阻力系数几乎不再减小,经过48轮迭代达到优化算法收敛条件,优化停止,结果如图7所示。优化算例1与初始构型俯仰力矩配平结果的力系数对比见表2。经过优化之后,阻力系数从初始的0.034 332下降至0.033 091,阻力减小3.61%,并且优化最终构型实现了纵向俯仰力矩的配平。迎角从初始的1.866°增加到了2.000°(迎角设计变量的边界)。图8给出了优化算例1与初始构型两者的表面压力分布的对比,初始构型机翼上翼面的等压线呈现汇聚、挤压的形态,表明激波较强;经过优化后,等压线趋于平行分布,表明通过优化消除了激波。

表1 优化算例1

表2 初始构型和有无短舱设计变量优化结果对比

图7 优化算例阻力系数收敛情况

图8 优化算例1与初始构型的表面压力分布的对比

5.3 考虑短舱位置变化的单点优化设计

在优化算例1的基础上增加短舱安装位置作为设计变量,将该优化算例命名为优化算例2,短舱的设计变量为下沉量、前伸量以及内偏角。优化模型见表3。

表3 优化算例2

优化算例2经过42轮迭代阻力系数几乎不再变化,经过57轮迭代达到优化算法收敛条件,优化停止,如图7所示。由于设计变量增加,设计空间维度提高,寻优难度增加,与优化算例1相比,优化算例2收敛需要更多的迭代计算。从表2的优化结果可以看出,经过优化之后,阻力系数从初始的0.034 332下降至0.330 20,阻力减小3.82%,迎角从初始的1.866°增加到了2.000°(迎角设计变量的边界)。与优化算例1相比,优化算例2的阻力系数进一步减小0.000 07,在实际工程应用中,大型民机阻力系数每下降0.000 1,载质量可增加7.56%[18]。根据文献[19],在民用客机标模(CRM)多学科分析中,减小0.000 1的阻力系数,相当于结构质量减小842 kg,这说明考虑短舱安装参数的优化设计结果具备一定的实际工程意义。图9给出了优化算例2与初始构型两者的表面压力分布的对比,与优化算例1类似,经过优化之后减弱了机翼上翼面的激波。

图9 优化算例2与初始构型的表面压力分布对比

短舱设计变量变化为内偏角增加0.201 5°,前伸量减小0.210 8 m,下沉量减小0.2 m。全机展向升力系数分布的对比如图10所示,可以看出,经过优化之后,机身所在区间载荷均比初始构型高,主要是通过优化使得机身迎角增加,增加了机身载荷,在全机升力系数一定的情况下,减小机翼所承担的载荷,减小激波阻力,更容易得到无激波的设计。通过扭转角和机翼几何型面的变化来改变机翼的载荷,降低内翼段载荷,增加外翼段载荷,进一步减小诱导阻力。

图10 两个优化算例与初始构型的展向单位化升力系数分布对比

(a)10.5%半展长(A剖面)(b)16.8%半展长(B剖面)(c)29.5%半展长(C剖面)

(d)34.1%半展长(D剖面)(e)49.2%半展长(E剖面)(f)90.8%半展长(F剖面)图11 两个优化算例与初始构型的剖面压力系数对比

对比优化算例1和优化算例2可以看出,短舱内偏角增加会减小机翼内翼段载荷,增加外翼段载荷,有利于减小诱导阻力。两个优化算例与初始构型在机翼不同位置的剖面压力系数的对比如图11所示。可以看出,通过调整短舱安装参数,减小前伸量和下沉量,降低了机翼剖面上翼面压力分布的头部负压峰值,并使得短舱所在处机翼剖面的下翼面压力分布更加光滑,如图11c所示。上翼面压力分布头部负压峰值的降低,减小了当地流动速度,避免了由于流动急剧加速导致的激波阻力增加;下翼面压力峰值的降低带来了额外的升力增量(等效于压力分布包围的面积)。优化设计结果综合论证了优化算例2具有减小机翼短舱相互干扰的效果。