图1 基于K r i g i n g 模型的优化设计算法框架 基于K r i g i n g 模型的翼型多目标气动优化设计研究
任庆祝，宋文萍
(西北工业大学翼型、叶栅空气动力学国防科技重点实验室，陕西西安710072
)
摘 要：在翼型气动优化设计中引入K r i g i n g 代理模型，发展了一套高效、稳定的气动优化设计程序。采用拉丁超立方试验设计方法在设计空间内构造一系列样本点，通过求解二维可压缩的雷诺平均N S 方程(R A N S )得到其响应值来建立初始K r i g i n g 模型。优化设计采用H i c k s -h e n n e 函数对翼型几何外形进行参数化表示，以阻力极小化为设计目标，考虑面积、升力、力矩等约束条件，通过算例证明，发展的优化设计方法不仅可行，而且具有高效稳定的特性。与传统的优化设计方法比较，大大减少了设计时间。 关键词:K r i g i n g 模型；试验设计；翼型优化设计
中图分类号:V 211.3 文献标识码:A 文章编号:1671-654X (2009)03-0077-06
引言
翼型的选择和设计一直是飞行器气动设计中的最
重要的环节之一，翼型选择的合适与否直接决定了飞行器的气动品质是否优良。当前翼型设计方法主要有数值优化设计
[1]
、反设计等等，传统的翼型优化设计耗
时多，主要体现在对目标函数的评估时需要大量调用流场求解程序，近些年来，应用代理模型技术
[5]
成为许
多气动优化设计者的首选，一些算例表明基于代理模型的设计方法不仅高效而且可靠，已逐渐在工程上得
到应用。
一般来说，应用于近似技术的代理模型有二次响应面模型,K r i g i n g 模型，径向基函数模型，人工神经网络模型等等，基于各种代理模型的理论研究和数值方法验证在许多期刊、专著上都有涉及
[2,4,5]
，本文采用
了精度较高、使用较为广泛的K r i g i n g 模型[6]
,K r i g i n g
模型起源于地理空间统计学，是一种估计方差最小的无偏估计模型，通过相关函数的作用，具有局部估计的特点，它可以较好的预估未知点处函数值的分布情况，进而可以代替原有的目标函数分析模型，在翼型气动优化设计上具有很大的应用价值。
1 基于K r i g i n g 模型的优化设计算法框架
基于代理模型的设计优化方法本质上是通过一些试验设计方法，分析建立响应量关于设计变量的近似模型，通过这些近似模型来研究最优设计问题。其步
骤为：首先对优化问题进行合理的数学描述，通过几何参数化确定设计变量及其取值范围，确定设计目标；然后采用试验设计方法在设计空间给出一定数量的初始样本点分布；根据初始样本点信息构造K r i g i n g 模型；对代理模型进行优化分析，寻最佳的设计变量取值；对搜索到的局部最优值，将它加入到初始样本点中，重新构造模型，如此迭代计算，直到收敛到全局的最优点。图1是基于K r i g i n g 模型的优化设计算法框架。
2 试验设计与近似模型
收稿日期:2008-04-10 修订日期:2008-07-29
基金项目：国家自然基金项目(90605004)；国家863高技术研究发展计划(2007A A 05Z 448) 作者简介：任庆祝(1984-)，男，陕西西安人，硕士研究生，研究方向为理论与计算流体力学。
第39卷 第3期航空计算技术
V o l .39N o .32009年5月
A e r o n a u t i c a l C o m p u t i n g T e c h n i q u e
M a y .2009
2.1 试验设计
为了能充分反映设计空间的特性，构造模型之前需选择一些有代表性的样本点，而样本点一般依据试
验设计[3]
的方法选择。试验设计方法是有关如何科学合理的安排试验的数学方法，常用的实验设计方法有全析因设计、正交设计、中心复合设计、均匀设计、拉丁方设计等，本文中采用使用较为广泛的拉丁超立方设计
[4,6]
。
2.2 构造K r i g i n g 模型
给定一组m 满足标准化条件的样本点X=[x 1
,
x 2
,…,x m ]T ，其中x i
是n 维行向量，其维数n 为设计
变量的个数，对应的响应为Y =[y 1,y 2,…,y m ]T
，其中
y i 是q 维行向量。第i 维响应y i
，
可以描述为:y i =f (βi ,X )+z i
(X )(1)
式中f (βi ,X )为回归模型部分(通常是多项式)，它是一个确定性部分,βi 为回归系数列向量;f (βi ,X )=[f 1(X ),…,f p
(X )]βi =f (X )βi (2)z i
(X )为一随机过程(i =1,2,…,q )有如下统计特性:E [z i (x )]=0(3)v a r [z i (z )]=σ2
i
(4)E [z i (x k
),z i (x j
)]=σ2
i R
(θ,x k
,x j
)(5)
式中R (θ,x k
完美前传
,x j )为两个样本点之间的相关函数，它是由用户定义的，在本文中取为G a u s s 函数;
R (θ,x k ,x j )=e x p [-∑n
k =1
θk x i k -x j
k 2
]
(6)
式中n 为设计变量的个数,x i k
和x j
k 为样本点x i
和x j
的第k 维分量;θk 为未知参数，在不同方向可以取不同值。根据模型要求预测模型的方差最小，可以得出待
测点的响应估计值为:
y (x )=f T β＊+r T R -1(R-F β＊)=f T β＊+r T γ＊(7)式中:
R i j =
R (θ,x i
,x j
),i ,j =1,2,…,m ,r (x )=[R (θ,x 1
,x ),…,R (θ,x m
,x )]
T
其中,x 为待测点,β＊
,σ2
可根据如下表达式估计;
β＊
=(F T
R -1
F )-1
F T
R -1
Y (8)σ 2
=
(Y -F β＊
)T
R -1
(Y -F β＊
)m
(9)
可以看出，在样本点一定的情况下,β＊
,γ＊
的值是
固定不变的，所以求解待测点的响应估计值时只要计算f (x ),r (x )就可以了。
对于θ的值的确定可以用如下方法，极大化该似然函数:
L n (β＊,σ2,θ)=-[m l n (σ 2
)+l n R ]/2
(10)
其中,θ>0，最好的模型即求解该无约束非线性最
优问题得到的。由K r i g i n g 模型预估未知点处的响应
值，其偏差或不确定性可通过如下的表达式得出，均方误差(M S E ):
s 2=φ(x )=F [y (x )-y (x )2
]=
σ2(1+u T (F T R -1F )-1u -r T R -1
城市雕塑网
r (11)u =F T
R -1
R-f
(12)
图2给出了K r i g i n g 初始模型近似一个一般函数的
示意图，可以看出在x =0.8附近由于没有样本点经过，预估值与函数真实值误差较大，这由最下方的均方根误差线也可以看出。因此，如图3所示，优化过程中需要在这个地方增加新的样本点以提高模型的预估精度。
图2 初始K r i g i n g 模型预估示意图
图3 更新后的K r i g i n g 模型预估示意图
3 模型可行性验证
K r i g i n g 模型建立后，要另取样本点来验证模型的精度，以保证模型的有效性，经过验证后的称K r i g i n g 模型才可以用来代替分析模型进行近似分析。一般可
利用均方根误差s (R M S E )和最大误差来验证K r i g i n g 模型。
均方根误差上文已经提到，最大误差为：在k 测试样本点中寻相对误差最大的一个值来检验模型的预估精度。
m a x e r r o r = y i -y i
y i
(i =1,…,k )(13)
本文采用了一种称为“标准化的交叉验证残值(s t a n d a r d i z e d c r o s s -v a l i d a t e d r e s i d u a l )”的方法
[6]
，先假
定模型有99.7%的置信水平，则相应的标准化的交叉
验证残值应位于[-3,+3]区间内，从而保证建立的
·78· 航空计算技术 第39卷第3期
模型具有足够的拟合精度:
y (x i )-y -i (x i
)s -i (x i
)(14)
4 收敛标准
真实目标函数的最小值下降的相对很小:(f m i n (x )i -f m i n (x )i -1)
f m i n (x )i -1
≤ε(15)
5 算例
5.1 测试函数
F =2x 2
1+3x 2
2-0.8s i n (2πx 1)-1.2c o s (3πx 2) x 1,x 2 [-1.024,1.024]F m i n =-1.889084,(x 1=0.2218,x 2=
0.0000)(16)
表1 测试函数
劳动保险条例实施细则修正草案
模型更新最小值点
K r i g i n g 预估
真实函数值
预估误差(%)0x 1=
0.2421-1.89878-1.880320.982x 2=-0.00501x 1=0.2260-1.88484-1.888370.187x 2=-0.00262
x 1=0.2209-1.88926-1.88907
0.100
x 2
=0.
0002图4 真实的函数响应值分布
本算例采用一个具有多个峰值的解析函数，其解析最小值为y m i n =-1.889084，在x (0.2218,0.0000)处，构造初始模型采用拉丁超立方设计使用51个样本点，图4、5给出了用K r i g i n g 模型近似真实函数的情
况，由图中可以看出模型对真实函数的预估总体上是较好的，从图中也可以看出该模型对全局极小值点的预估比较完全，通过表1也可以看出经过两次模型更新后，得到全局的最小值(y =-1.88907)与解析最小
值非常接近
。
图5 预估的函数响应值分布
5.2 翼型气动优化设计
本文采用H i c k s -H e n n e 形状函数来描述翼型的几何扰动，新的翼型外形为基本翼型加扰动构成，其表示
如下:
y=y b a s e +∑n v
j =1
δi f i
(x )(17)其中,10个H i c k s -H e n n e 函数是:f 1(
x )=x 0.25
(1-x )e -20x
(18)f i
(x )=s i n 3
[πx e (i )
] (i =2,3,4,5)(19)e (i )=l n (0.5)l n (x i
) (i =2,3,4,5)(20)x i =
0.2,0.4,0.6,0.8 (i =2,3,4,5)(21)
office小f i +5(x )=-f i (x )(22)式(17)中,n v 为设计变量的总数，形函数f (x )前5个分布于翼型上表面，后5个分布于翼型下表面。式中各形函数的系数δ即为设计变量，其取值范围可在[-4.＊t h i c k 0/1.e +2,4.＊t h i c k 0/1.e +2]；本文采用复合形法作为数值优化方法，构造初始K r i
g i n g 模型采用了121个样本点(每个样本点响应值采用R A N S 方程求解)，表2，表3给出了优化设计的结果。
单目标的翼型减阻设计算例:
算例1
本文将代理模型技术应用于翼型气动优化设计，以面积、升力、力矩等作为约束条件，进行减阻设计，基本翼型为R A E 2822，迎角2.3°，马赫数0.73，雷诺数6.5E +6，目标函数、约束条件分别表示为:
D e s i g n c o n d i t i o n :M a =
0.73,α=2.3°,R e =6.5E +6M i n i m i z e C d
A r e a ≥A r e a 0＊0.99S u b j e c t t o c l ≥C 10＊0.99
1.1·
79·2009年5月任庆祝等：基于K r i g i n g 模型的翼型多目标气动优化设计研究
表2 第一次优化后的结果
参数基本翼型优化翼型(K r i g i n g 预估)优化翼型(R A N S 方程)变化量(%)误差(%)均方根误差(%)
C l 0.767560.769740.76548-0.27080.55261.6333C d 1.4336E-21.0343E-21.1315E-2-21.0779.39920.0906C m -1.0221E-1-9.8628E-2-9.9280E -22.86730.66040.5300A 7.7723E-2
7.7160E-2
-0.7244
表3 第二次优化的结果
参数基本翼型
优化翼型(K r i g i n g 预估)优化翼型(R A N S 方程)变化量(%)误差(%)均方根误差(%)
C l
0.76756
0.763930.76393-0.47371.29101.3213
C d 1.
4336E-21.1424E-21.1284E-2
-21.290
1.2269
0.0742
C m -1
.0221E-1-1.0256E-1
-9.7809E -24.30584.6375
0.4297
A 7.7723E-2
7.7225E-2
-0.
6408
图6 模型的交叉验证残值C L 图7 模型的交叉验证残值C D 图8
优化前后的翼型几何外形比较
图9 优化前后的压力分布比较 图10 第一设计点优化前后的表面压力 图11 第二设计点优化前后的表面压力
算例1构造初始K r i g i n g 模型安排了121个试验
点，调用了流场求解程序121次，图6、7分别为预估响应值为升力系数和阻力系数时模型的交叉验证残值，显然该模型的预估精度是满足使用要求的。表2,3给出了分别在第一次、第二次优化后结果，第一次优化后阻力系数下降了21.07%，力矩系数绝对值下降了2.86%，升力系数略有下降，第二次优化后阻力系数下降了21.29%，力矩系数绝对值下降了4.30%，升力系数也下降0.47%，但仍在约束范围内。同时从表中也能看出模型对最优点的预估精度也提高了，从图8、9可以看出优化后的翼型(O P T F O I L 1,O P T F O I L 2)激波强
度与原始翼型相比明显减弱，翼型的前缘几何外形有
尺度函数了明显变化。
多目标的翼型减阻设计
对于多目标的优化问题，本文采用“统一目标函数化”[1]
进行了处理，使其转化为一个单目标的优化问题；假设第一设计点和第二设计点的目标函数分别为F 1(x ),F M 2(x )，分别给其加权因子ω1,ω2，例如(0.35,0.65)，基本翼型在两个设计状态下的阻力分别为C d 01,C d 02
令
α1=0.8C d 01,α2=0.8C d 02
β1=1.2C d 01,β2=1.2C d 02
(21)
·80· 航空计算技术 第39卷第3期
则规格化后的目标函数为:F ′i
(x )=F i (x )-αi
βi -αi (i =1,…,q )(22)
“统一的目标函数”为
F (x )=∑q
i =1
ωi F ′i
(x )(23)
算例2
D e s i g n c o n d i t i o n (1):M a =0.73,α=1.3420
,R e =6.5E +6
D e s i g n c o n d i t i o n (2):M a =
0.73,α=1.8420
,R e =6.5E +6
M i n i m i z e C d 1,C d 2
A r e a ≥A r e a 0＊0.99S u b j e c t t o C l 2≥C l 02＊0.
99C m 2
C m 02
-1.01≤0表4 第一设计点优化结果
参数基本翼型优化翼型(K r i g i n g 预估)优化翼型(R A N S 方程)变化量(%)误差(%)均方根误差(%)
C l 0.583180.576080.56922-2.39121.19211.1323C d 9.5199E-39.0510E-39.2501E-3-2.83462.19920.0271C m -9.9627E-2-9.8712E-2-9.7901E -21.73290.82220.3289A 7.7723E-2
7.7034E-2
-0.8868
表5 第二设计点优化结果
参数基本翼型优化翼型(K r i g i n g 预估)优化翼型(R A N S 方程)变化量(%)误差(%)均方根误差(%)
C l
0.68112
0.67800
0.66962
-1.6889
1.2358
1.1182
C d
1.1269E-2
9.7545E-31.0031E-2-10.9812.83720.0402C m -1
.0078E-1-9.7865E-2
-9.6943E -2
3.8102
0.9422
0.3340
算例2同样采用了242个初始样本点构造K r i g i n g 模型，表4,5给出了优化结果，经过一次优化后两个设计
状态的阻力系数分别下降了2.83%,10.98%，同时升力、力矩、面积都有所下降，与参考文献[1]比较，考虑到模型近似的误差影响，这个优化结果也是可接受的。
6 结论
贝司大师本文将K r i g i n g 代理模型技术与数值优化方法结合，进行了应用于翼型气动优化设计的研究，通过算例证明，这套方法不仅可行而且高效、实用，具有较强的工程应用价值，值得注意的是，约束条件、试验设计的样本点分布的正交性、均匀性特点，以及数值优化方法的全局性和稳定性也对设计结果具有重要影响，下一步，可针对这些问题加以改进，针对复杂外形气动优化设计作深入研究。参考文献:
[1] 刘战合，宋文萍，宋笔锋.基于复合形方法的翼型优化设
计研究[J ].西北工业大学学报,2004,22(6):795-799.[2] 熊俊涛，乔志德，韩忠华.基于N a v i e r -S t o k e s 方程跨声速
翼型和机翼气动优化设计[J ].空气动力学学报,2007,25
(1):29-33.
[3] A n t h o n y AG i u n t a ,S t e v e n FWo j t k i e w i c z J r ＆M i c h a e l SE l -d r e d ,O v e r v i e w O f M o d e r nD e s i g nO f E x p e r i m e n t s M e t h o d s F o r C o m p u t a t i o n a l S i m u l a t i o n s [A ].A I A A 2003-0649.[4] S h i n k y uJ e o n g ,K a z u o m i Y a m a m o t oa n dS h i g e r uO b a y a s h i ,
K r i g i n g -B a s e dP r o b a b i l i s t i cM e t h o df o rC o n s t r a i n e dM u l t i -O b j e c t i v eO p t i m i z a t i o nP r o b l e m[A ].A I A A p a p e r 2004-6437.
[5] 张科施，韩忠华，李为吉.基于近似技术的高亚声速运输
机机翼气动/结构优化设计[J ].航空学报,2006,27(5):810-815.
[6] S r e nNL o p h a v e n ,H a n s B r u u n N i e l s e n ,J a c o bSn d e r g a a r d ,
D A C
E :aM a t l a b K r i g i n gT o o l b o x [A ].T e c h n i c a lR e p o r t I M M -T R-2002-12.
[7] D O N A L D R J O N E S ,M A T T H I A SS C H O N L A U ,WI L L I A M J
WE L C H .E f f i c i e n t G l o b a l O p t i m i z a t i o no f E x p e n s i v eB l a c k -B o x F u n c t i o n s [J ].J o u r n a l o f G l o b a l O p t i m i z a t i o n 13:455-492,1998.
[8] M a s a t oS e k i s h i r o ,G e r h a r dV e n t e r a n dV l a d i m i r B a l a b a n o v .
C o m b i n e dK r i g i n g a n dG r a d i e n t -B a s e dO p t i m i z a t i o nM e t h o d [A ].A I A A 2006-7091.
·
81·2009年5月任庆祝等：基于K r i g i n g 模型的翼型多目标气动优化设计研究
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议