基于Kriging方法的整流罩气动外形设计优化①

杨希祥，周 张，彭 科

(1.国防科技大学航天科学与工程学院，长沙 410073;2.北京理工大学 宇航学院，北京 100081)

0 引言

气动外形设计在运载火箭总体设计中占有极其重要的地位。运载火箭外形一般由头部、箭身和尾舵(小型运载火箭)等部分组成。头部外形设计即整流罩外形设计对飞行性能影响至关重要，是外形设计的关键［1］。运载火箭总体设计需重点考虑的气动阻力、表面传热特性和有效载荷容积等因素，都与整流罩外形息息相关。整流罩外形设计在国外飞行器设计中历来受到高度重视［2－4］，国内却鲜有文献对这一问题进行研究。

气动外形设计优化过程往往需要多次调用气动参数计算模型，直接采用高精度计算模型将面临严重计算复杂性问题，近似建模方法应运而生。近似建模方法基本思想是用一个简单逼近函数近似替代复杂高精度实际性能分析模型。根据近似函数所能模拟设计空间的大小，近似建模方法可分为局部近似方法、全局近似方法和中等范围近似方法。全局近似方法中的Kriging函数法以具有样本点处无偏估计、良好的高阶非线性拟合能力和灵活的近似模型参数选择等优点，在近似建模领域得到广泛应用。Meunier和Laurenceau 等采用其建立机翼气动力近似模型［5－6］;Han等提出一种基于Kriging函数的变复杂度近似建模方法，并应用于气动力近似建模［7］;Ahmed等分别研究了Kriging函数在高超声速钝锥体气动力和气动热近似建模中的应用［8］。

本文研究小型运载火箭整流罩气动外形设计优化问题，近似模型构造选用Kriging函数法。

1 Kriging函数法及其仿真实现

1.1 参数估计与预测

假定n维输入变量X与输出响应y之间的函数关系y(X)未知，给定组s训练样本数据(X，y)=［(X1，y1)，…，(Xs，ys)］T，结合线性回归参数模型和随机过程非参数模型，采用Kriging函数可构造y(X)的近似模型 y^(X)如下［9］:

式中 fT(X)=［f1(X)，…，fe(X)］为回归函数;β=［β1，…，βe］T为回归系数;z(X)为随机过程。

工程应用中，线性回归参数模型常采用二阶多项式回归函数，

随机过程非参数模型考虑为真实响应与假设线性模型的偏差，常采用均值为零、高斯相关函数的随机过程模型，其协方差函数为

随机变量 y=y(X)和 Ys=［y(X1)，…，y(Xs)］T服从如下多元正态分布:

其中:

式中 fT为e维回归函数行向量;rT为s维相关函数行向量;F为s×e维回归函数矩阵;R为s×s维相关函数矩阵。

y Ys服从如下条件分布:给定训练样本数据Ys下，y的最小均方误差估计为yYs的条件均值:

要得到式(6)的Kriging函数近似模型，需基于训练样本数据估计近似模型参数β、θ，通常采用两步估计法，β与θ相互独立，即满足下式:

第一步估计是假定已知相关函数模型及参数θ，得到β的广义最小二乘估计(θ):

第二步估计是得到θ的某个估计θ^，然后将θ^代替θ，得到β^(θ^)作为β的估计。相关参数θ的估计方法有交互验证估计方法和极大似然估计方法，训练样本点数较少时，可选择交互验证估计方法，否则选择极大似然估计方法。

相关参数θ的交互验证估计值θ^为所有样本训练点处的均方预估误差的最小值点:

采用极大似然估计方法时，假设训练样本Ys=［y(X1)，…，y(Xs)］T服从多元正态分布:

取似然函数FMLE为训练样本数据Ys概率密度f(Ys):

不考虑常数项的对数似然函数为

将式(15)代入式(14)，得到只与相关参数θ相关的对数似然函数:

式(18)所示的近似模型不仅可以预估输入变量X和输出响应关系y(X)，而且还可以预估输入变量X处一阶偏导数:

1．2 仿真实现

实际应用中，由于训练样本数量大，致使高维相关函数矩阵的求逆运算计算复杂，同时训练样本点选择不合理会导致相关函数矩阵病态，无法计算逆矩阵。由于相关函数矩阵对称，属稀疏矩阵，为避免直接计算产生较大计算误差，对相关函数矩阵进行Cholesky分解;为避免不同的设计变量量级不同产生计算误差，对训练样本数据进行标准化处理，具体步骤如下:

(1)对训练样本数据进行标准化处理，得到标准空间［－1，1］s上的训练样本数据

(4)对回归函数矩阵F进行正交转化并做QR分解:QG=C－1F。

(5)计算基于标准化训练样本数据的Kriging函数近似模型:

(6)计算式(11)均方预估误差值fCV(θ)或式(17)对数似然函数值LMLE(θ)。

(7)更新相关参数比例参数θ，重复执行步骤(3)～(6)，直到找到目标函数极小值。

(8)基于训练样本数据的Kriging近似模型:

2 整流罩气动外形设计优化问题

2．1 设计变量选取

研究的运载火箭整流罩前罩基线外形如图1所示。为灵活表示整流罩气动外形，便于选取控制变量，采用易修改性和局部调整特性好的非均匀有理B样条曲线(Non-Uniform Rational B-Spline，NURBS)对整流罩外形进行参数化表示。

控制点选取如图2所示，除首末2点外，另取4个控制点，共6个控制点。其中，控制点2和控制点1具有相同横坐标(均为0)，目的是保持NURBS曲线光滑性，控制点3距离头部顶点很近，目的是增强外形多样性，控制点5与控制点6具有相同纵坐标，目的是保持后端光滑性［4］。根据NURBS曲线性质，改变控制点坐标即可方便调整整流罩外形，为此，选取控制点坐标为设计变量，共4个，具体如下:控制点3的横、纵坐标分别为x1和x2，控制点4的纵坐标为x3，控制点5的横坐标为x4，即控制变量:

图1 整流罩基线外形Fig.1 Primary shape of fairing

图2 控制顶点选取Fig.2 Selection of control points

设计变量边界值的确定对优化设计速度和结果有重要影响。为加快收敛速度，同时确保得到全局最优解，首先采用DATCOM气动特性工程估算程序与Fay和Riddell热流密度计算公式，结合Powell法，求得粗略最优解，用于收缩设计空间，然后采用基于Kriging函数近似模型的方法进一步精确求解。

DATCOM是美国空军实验室研究的一套用于导弹气动力工程估算的程序，对整流罩这样的标准细长体预估误差在10%～15%，对DATCOM源程序进行校核和修正，形成核心计算程序，编制输入输出文件，嵌入优化流程对迭代过程气动力进行计算。Fay和Riddell热流近似计算公式用于求解轴对称体驻点区热流密度:

式中 下标w、e、s分别表示壁面条件、边界层外缘条件、驻点条件;qws为驻点热流密度;ρw、ρs为密度;μw、μs为粘性系数;hs、hw为比焓;hd为离解比焓;RN为钝头半径。

为便于计算，Le和hd用下式近似计算［10］:

2.2 目标函数和约束条件

运载火箭大气层飞行段阻力对运载能力有重要影响，而气动阻力主要集中在头部。为此，以一级飞行段平均阻力最小为目标函数，选取基线方案3个典型弹道特征点，阻力系数最大点、动压最大点和一级发动机关机点，3点相关参数如表1所示。

表1 弹道特征点参数Table 1 Trajectory characteristic parameters

引入加权因子对目标函数进行处理，约束条件包括3个设计点驻点热流密度最大值、整流罩内部容积及设计变量约束:

式中 Di、qwi(i=1，2，3)分别表示设计点i阻力和驻点热流密度;qw，max为允许热流密度最大值;V表示整流罩内部容积;Vb表示基线方案整流罩容积值，设计变量边界根据初步优化值和设计变量实际允许的边界值得到。

3 设计优化结果与分析

采用Kriging函数建立运载火箭整流罩气动性能计算近似模型，采用正交表型均匀拉丁超立方试验设计方法［11］选取样本点，样本点选取36个。

样本训练采用非结构网格和基于N-S方程的数值计算方法，湍流模型采用k-ε两方程模型，空间离散采用上风差分格式中的通量差分分裂(FDS)格式，时间分裂采用LU-SGS方法，由于3个特征设计点训练样本总数较多，采用并行计算产生气动力参数。优化方法选取Powell法，初值设定采用2.1节所述方法。

式(25)中，w1、w2和 w3分别取为 0.2、0.5、0.3。经4次迭代优化过程收敛，优化方案目标函数和约束条件如表2所示，表中同时给出了基线方案和在3个特征点分别以阻力最小为目标函数设计结果对应的目标函数和约束条件值。

优化前后整流罩外形(含计算网格)对比如图3所示。

特征点2处，优化方案和基线方案表面压力和热流分布对比如图4所示。

表2 整流罩外形优化设计结果Table 2 Optimization design results of fairing shape

图3 基线方案与优化方案外形对比(含计算网格)Fig.3 Comparison of shape between basic scheme and optimal scheme(computation grid)

图4 基线方案和优化方案轴向参数分布对比Fig.4 Comparison of axial parameters between basic scheme and optimal scheme

由图4可见，由于头部钝度减小，优化方案在头部具有更小表面压力，但在距头部约45～650 mm处，优化方案表面压力大于基线方案;优化方案母线光滑过渡，表面压力分布脉动较小，有利于整流罩结构设计;由图4还可看出，头部钝度减小使得优化方案驻点热流密度较高，但满足设计约束。

4 结论

(1)采用Kriging函数建立整流罩气动参数计算近似模型，可有效避免整流罩气动外形设计优化计算复杂性问题，在保证优化精度的前提下，实现较高的优化效率。

(2)采用建立的整流罩气动参数计算近似模型，实现了运载火箭整流罩气动外形设计优化问题，优化方案平均阻力比基线方案减小22.2%，各项约束条件得到良好满足。

(3)研究成果对于方案论证与方案设计阶段飞行器气动外形设计优化研究具有重要借鉴意义。

［1］龙乐豪.总体设计(上)［M］.北京:宇航出版社.1989.

［2］Lee Young Ki，Lee Jae-Woo，Byun Yung-Hwan.The design of space launch vehicle using numerical optimization and inverse method［C］//17th AIAA Applied Aerodynamics Conference.Norfolk，VA，1998:757-767.

［3］Deeppark N R，Ray T，Boyce R R.Evolutionary algorithm shape optimization of a hypersonic flight experiment nose cone［J］.Journal of Spacecraft and Rockets，2008，45(3):428-436.

［4］Lee Jae-Woo，Min Byung-Young，Byun Yung-Hwan，et al.Multipoint nose shape optimization of space launcher using response surface method［J］.Journal of Spacecraft and Rockets，2006，43(1):137-146.

［5］Meunier M.Simulation and optimization of flow control strategies for novel high-lift configurations［J］.AIAA Journal，2009，47(5):1145-1157.

［6］Laurenceau J，Sagaut P.Buliding efficient reponse surfaces of aerodynamic functions with Kriging and Cokriging ［J］.AIAA Journal，2008，46(2):498-507.

［7］Han Z H，Ralf Z，Gortz S.A New Cokriging method for variable-fidelity surrogate modeling of aerodynamic data［C］//48th AIAA Aerospace Sciences Meeting.Orlando，Florida，2010:1-17.

［8］Ahmed M Y M，Qin N.Metamodels for aerothermodynamic design optimization of hypersonic spiked blunt bodies［C］//48th AIAA Aerospace Sciences Meeting.Orlando，Florida，2010:1-17.

［9］Matheron G.Principles of geostatistics［J］.Economic Geology，1963(58):1246-1266.

［10］王国雄.弹头技术(上)［M〛.北京:宇航出版社，1993.

［11］张润楚，马长兴.正交表型均匀LH设计和抽样［J］.应用概率统计，2001，17(3):243-254.