段焰輝 吳文華 范召林 羅佳奇

1)(中國(guó)空氣動(dòng)力研究與發(fā)展中心,計(jì)算空氣動(dòng)力學(xué)研究所,綿陽(yáng) 621000)2)(北京大學(xué)工學(xué)院,北京 100871)

基于本征正交分解的氣動(dòng)優(yōu)化設(shè)計(jì)外形數(shù)據(jù)挖掘?

段焰輝1)?吳文華1)范召林1)羅佳奇2)?

1)(中國(guó)空氣動(dòng)力研究與發(fā)展中心,計(jì)算空氣動(dòng)力學(xué)研究所,綿陽(yáng) 621000)2)(北京大學(xué)工學(xué)院,北京 100871)

(2017年7月2日收到;2017年7月18日收到修改稿)

氣動(dòng)外形的全局優(yōu)化設(shè)計(jì)會(huì)產(chǎn)生大量的過(guò)程數(shù)據(jù),其中隱含的設(shè)計(jì)知識(shí)具有較高的挖掘價(jià)值.數(shù)據(jù)挖掘有助于獲取直觀、可定性描述的設(shè)計(jì)知識(shí).本文采用基于本征正交分解的數(shù)據(jù)挖掘方法從氣動(dòng)優(yōu)化設(shè)計(jì)的過(guò)程數(shù)據(jù)中獲取設(shè)計(jì)知識(shí),數(shù)據(jù)挖掘?qū)ο鬄榭缫羲賶簹鈾C(jī)轉(zhuǎn)子葉片NASA Rotor 37的優(yōu)化過(guò)程數(shù)據(jù),該數(shù)據(jù)由基于粒子群方法的絕熱效率最大化優(yōu)化設(shè)計(jì)產(chǎn)生.結(jié)果表明:基于本文數(shù)據(jù)挖掘方法獲取的設(shè)計(jì)知識(shí)能夠直接反映氣動(dòng)外形的變化規(guī)律,為葉片的氣動(dòng)外形設(shè)計(jì)提供參考;數(shù)據(jù)挖掘的設(shè)計(jì)知識(shí)成功地驗(yàn)證了優(yōu)化設(shè)計(jì)結(jié)果的有效性.

數(shù)據(jù)挖掘,本征正交分解,氣動(dòng)優(yōu)化設(shè)計(jì),跨音速

1 引 言

氣動(dòng)外形優(yōu)化設(shè)計(jì)是航空航天領(lǐng)域的重要學(xué)科之一,相較傳統(tǒng)的設(shè)計(jì)方法,具有成本低、周期短的突出優(yōu)點(diǎn).隨著氣動(dòng)外形優(yōu)化設(shè)計(jì)技術(shù)的發(fā)展,隨機(jī)搜索類方法因?yàn)槟軌蚩朔鹘y(tǒng)優(yōu)化方法局部性的缺陷,得到了日益廣泛的應(yīng)用.這類優(yōu)化方法在優(yōu)化設(shè)計(jì)中將產(chǎn)生大量的過(guò)程數(shù)據(jù),但是只有最終的優(yōu)化結(jié)果被保留應(yīng)用.耗費(fèi)大量計(jì)算資源的過(guò)程數(shù)據(jù)并非毫無(wú)用處,其中隱含了與優(yōu)化相關(guān)的設(shè)計(jì)知識(shí),即設(shè)計(jì)變量與目標(biāo)函數(shù)的關(guān)系.若能獲取這些數(shù)據(jù)中的設(shè)計(jì)知識(shí)并加以利用,有助于揭示設(shè)計(jì)空間內(nèi)設(shè)計(jì)變量變化對(duì)目標(biāo)函數(shù)影響的一般規(guī)律,更深層次地分析認(rèn)識(shí)設(shè)計(jì)變量對(duì)目標(biāo)函數(shù)的影響機(jī)制,為類似的設(shè)計(jì)優(yōu)化積累設(shè)計(jì)信息和經(jīng)驗(yàn),對(duì)提高類似設(shè)計(jì)優(yōu)化的效率與可靠性具有重要的指導(dǎo)意義.本質(zhì)上,數(shù)據(jù)挖掘就是這樣一種從大量數(shù)據(jù)中獲取隱含知識(shí)的方法.本文將利用數(shù)據(jù)挖掘技術(shù)獲取優(yōu)化設(shè)計(jì)過(guò)程數(shù)據(jù)中的設(shè)計(jì)知識(shí).

目前得到應(yīng)用的數(shù)據(jù)挖掘方法較多,應(yīng)用領(lǐng)域極廣.在優(yōu)化設(shè)計(jì)方面,數(shù)據(jù)挖掘技術(shù)多用于多學(xué)科、多目標(biāo)優(yōu)化設(shè)計(jì)的設(shè)計(jì)知識(shí)挖掘[1].Jeong等[2]、Chiba和Obayashi[3]采用總變差分析和自組織映射方法對(duì)兩級(jí)入軌飛行器可飛回第一級(jí)的多學(xué)科、多目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果進(jìn)行了分析,實(shí)現(xiàn)了從非劣解集中選擇最終結(jié)果;Oyama等[4,5]采用基于本征正交分解(proper orthogonal decomposition,POD)[6]的數(shù)據(jù)挖掘技術(shù)對(duì)跨音速翼型的多目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果進(jìn)行了分析,將優(yōu)化結(jié)果分成了三類,并獲取了每類結(jié)果的設(shè)計(jì)知識(shí).郭振東等[7]對(duì)數(shù)據(jù)挖掘的應(yīng)用進(jìn)行了拓展,采用顯著變量識(shí)別、總變差分析及自組織映射等數(shù)據(jù)挖掘技術(shù)對(duì)設(shè)計(jì)空間進(jìn)行知識(shí)挖掘,得到的設(shè)計(jì)知識(shí)與優(yōu)化設(shè)計(jì)結(jié)果進(jìn)行了相互驗(yàn)證.汪偉等[8]將優(yōu)化設(shè)計(jì)的過(guò)程數(shù)據(jù)作為數(shù)據(jù)挖掘?qū)ο?采用K-Means聚類分析、粗糙集屬性重要性算法和決策樹(shù)方法對(duì)優(yōu)化過(guò)程數(shù)據(jù)進(jìn)行處理和分析,得到了葉片優(yōu)化設(shè)計(jì)變量與目標(biāo)函數(shù)的隱含規(guī)則.綜上可知,數(shù)據(jù)挖掘在氣動(dòng)外形優(yōu)化設(shè)計(jì)中的應(yīng)用尚處于探索階段,分析方法和對(duì)象呈現(xiàn)多樣化,具有重要的研究意義.

本文基于POD的數(shù)據(jù)挖掘技術(shù)應(yīng)用于優(yōu)化設(shè)計(jì)過(guò)程數(shù)據(jù)的處理,以獲取設(shè)計(jì)知識(shí).為了提高數(shù)據(jù)挖掘?qū)ο蟮目尚哦?也即數(shù)值精度,本文的優(yōu)化過(guò)程數(shù)據(jù)都由計(jì)算流體力學(xué)(computational fluid dynamics,CFD)方法確定.為了縮短優(yōu)化周期,采用同步并行粒子群方法(parallel synchronous particle swarm optimization,PSPSO)在大規(guī)模集群上進(jìn)行優(yōu)化.此外,當(dāng)數(shù)據(jù)挖掘?qū)ο鬄闅鈩?dòng)優(yōu)化設(shè)計(jì)的設(shè)計(jì)變量,結(jié)果的直觀性普遍較差,因?yàn)槎鄶?shù)參數(shù)化方法提取的設(shè)計(jì)變量都無(wú)法直接反映氣動(dòng)外形的特征.基于POD的數(shù)據(jù)挖掘可以直接對(duì)外形數(shù)據(jù)進(jìn)行分析,數(shù)據(jù)挖掘的結(jié)果能直接反映外形的變化,提供直觀的設(shè)計(jì)知識(shí),因而本文采用POD進(jìn)行數(shù)據(jù)挖掘,改善結(jié)果的直觀性.

本文首先介紹POD的基本原理及基于POD的數(shù)據(jù)挖掘方法.然后介紹本文采用的PSPSO氣動(dòng)外形優(yōu)化設(shè)計(jì)方法,以及三維葉片的參數(shù)化方法和網(wǎng)格變形方法;并采用PSPSO方法對(duì)跨音速壓氣機(jī)轉(zhuǎn)子葉片NASA Rotor 37進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì),提取過(guò)程數(shù)據(jù).最后采用POD對(duì)過(guò)程數(shù)據(jù)進(jìn)行挖掘,分析所隱含的設(shè)計(jì)知識(shí),并驗(yàn)證本文所提出的數(shù)據(jù)挖掘技術(shù)的有效性.

2 POD數(shù)據(jù)挖掘方法

數(shù)據(jù)挖掘方法一般可以分為三步:數(shù)據(jù)制造、數(shù)據(jù)預(yù)處理和數(shù)據(jù)挖掘.本文對(duì)優(yōu)化設(shè)計(jì)的過(guò)程數(shù)據(jù)進(jìn)行數(shù)據(jù)挖掘,數(shù)據(jù)制造對(duì)應(yīng)優(yōu)化設(shè)計(jì),數(shù)據(jù)預(yù)處理對(duì)應(yīng)制造數(shù)據(jù)的篩選和處理,數(shù)據(jù)挖掘?qū)?yīng)優(yōu)化氣動(dòng)外形的特征分析.本文的數(shù)據(jù)預(yù)處理主要是按目標(biāo)函數(shù)的優(yōu)劣對(duì)優(yōu)化過(guò)程數(shù)據(jù)進(jìn)行篩選,具體做法是:與原始外形的氣動(dòng)性能相比,將氣動(dòng)性能提高的過(guò)程數(shù)據(jù)篩選出來(lái),剔除氣動(dòng)性能降低的過(guò)程數(shù)據(jù).值得注意的是:優(yōu)化設(shè)計(jì)中得到的過(guò)程數(shù)據(jù)包含氣動(dòng)外形、流場(chǎng)等信息,在數(shù)據(jù)挖掘時(shí)并非所有數(shù)據(jù)都要使用,因而在預(yù)處理階段需要對(duì)這些數(shù)據(jù)予以區(qū)分.本文對(duì)氣動(dòng)外形進(jìn)行數(shù)據(jù)挖掘,因此在預(yù)處理階段只保留氣動(dòng)外形信息.

2.1 POD基本原理

原始POD方法[6]需要整個(gè)定義域中的解集,維度較高,對(duì)特征值求解的精度和穩(wěn)定性提出了較大挑戰(zhàn),難以在工程中得到應(yīng)用.目前多采用由Sirvoich提出的Snapshot POD方法[9],該方法用一組定義域中的快照代替整個(gè)定義域中的解集,提高了POD方法的穩(wěn)定性,降低了計(jì)算量,目前在氣動(dòng)研究中得到了大量應(yīng)用[10?15].下文中,若無(wú)特殊說(shuō)明,所涉及的“POD方法”都為“Snapshot POD方法”.

假設(shè)在某一定義域內(nèi)有一組快照:s1,s2,···,sm,其POD基為:φ1,φ2,···,φm,m為快照數(shù).POD基的求解方法如下.

首先求解特征值問(wèn)題,

其中R為快照集合的自相關(guān)矩陣,λ為特征值,a為特征向量.R定義為

由奇異值分解(singular value decomposition,SVD)確定R的特征值和特征向量后,POD基可由下式計(jì)算:

其中:φk表示第k個(gè)基,表示第k個(gè)特征值對(duì)應(yīng)的特征向量的第i個(gè)元素.

2.2 基于POD的數(shù)據(jù)挖掘

本文基于POD的數(shù)據(jù)挖掘與應(yīng)用主要采用文獻(xiàn)[12,13]介紹的基于POD的參數(shù)化方法,在此基礎(chǔ)上加入一些能夠反映物理特征的量來(lái)挖掘設(shè)計(jì)知識(shí).本文POD分析對(duì)象為優(yōu)化過(guò)程中的氣動(dòng)外形.

假設(shè)被分析的氣動(dòng)外形由表面網(wǎng)格點(diǎn)組成,定義一組向量g1,g2,···,gm,對(duì)于優(yōu)化設(shè)計(jì)的過(guò)程數(shù)據(jù)而言,這組向量應(yīng)該具有某種共性,比如優(yōu)化的目標(biāo)函數(shù)都高于初始值.每個(gè)向量中的數(shù)據(jù)為氣動(dòng)外形表面網(wǎng)格點(diǎn)三維坐標(biāo)的順序排列,如gi={X,Y,Z}T,X,Y,Z依次為網(wǎng)格點(diǎn)的三個(gè)坐標(biāo)值,并且按照網(wǎng)格點(diǎn)編號(hào)排成序列.若氣動(dòng)外形表面網(wǎng)格點(diǎn)數(shù)為n,則每個(gè)向量的維數(shù)為3n.由氣動(dòng)外形表面網(wǎng)格點(diǎn)坐標(biāo)組成的向量可以確定POD所需的快照:

其中,g0可以選為網(wǎng)格點(diǎn)坐標(biāo)向量的平均值,也可以是優(yōu)化的初始外形.每個(gè)外形快照減去g0,便于分析所有氣動(dòng)外形的擾動(dòng)特征.由SVD確定POD基后,所有快照都可由這些基表示:

數(shù)據(jù)挖掘的設(shè)計(jì)知識(shí)就隱含在POD基和系數(shù)(也稱為基的權(quán)重)中.基能夠反映氣動(dòng)外形變化的區(qū)域及相對(duì)大小,結(jié)合其系數(shù)就能得到氣動(dòng)外形變化的絕對(duì)大小,從而給出最終的氣動(dòng)外形變化規(guī)律.但是基較多(與快照數(shù)目相同),系數(shù)也較多,一般只需要對(duì)其中能夠反映主要規(guī)律的基和相應(yīng)系數(shù)進(jìn)行分析即可.如何篩選基和系數(shù),是基于POD數(shù)據(jù)挖掘方法的一個(gè)重點(diǎn).

2.2.1 POD基及基系數(shù)的篩選

引入基的能量的概念[4]對(duì)基進(jìn)行篩選,第i個(gè)POD基的能量Energy定義為

其中φi,j表示第i個(gè)POD基的第j個(gè)元素.由(6)式可知:能量其實(shí)就是每個(gè)外形的網(wǎng)格點(diǎn)位移大小的總體度量.再定義能量比重:

由能量比重可以確定每個(gè)POD基所占的比重,從而確定哪些基是獲取設(shè)計(jì)知識(shí)的關(guān)鍵.

基系數(shù)的篩選要解決兩個(gè)問(wèn)題:分析哪些基的系數(shù),分析哪些快照的系數(shù).本文研究中,只對(duì)篩選得到的POD基的系數(shù)進(jìn)行分析.快照則根據(jù)優(yōu)化目標(biāo)函數(shù)的大小,選擇靠前的部分快照,這樣可以保證得到的基的系數(shù)反映的都是最優(yōu)外形的設(shè)計(jì)規(guī)律.

2.2.2 設(shè)計(jì)知識(shí)的獲取方法

POD基及由POD基和基系數(shù)組合表示的氣動(dòng)外形,都是坐標(biāo)的擾動(dòng)信息,不易直接獲得設(shè)計(jì)知識(shí),需要將其變換成能夠反映物理特征的量.氣動(dòng)外形的變化規(guī)律實(shí)際上就是每個(gè)網(wǎng)格點(diǎn)位置的變化規(guī)律,這種變化不僅有大小信息,還有方向信息.確定這種變化量的最佳辦法是采用優(yōu)化時(shí)的參數(shù)化方法,根據(jù)參數(shù)化方法中每個(gè)網(wǎng)格點(diǎn)位置擾動(dòng)的施加方式(一般是在某個(gè)方向施加),反推網(wǎng)格點(diǎn)的位置變化.本文采用了前期研究中的參數(shù)化方法[16],擾動(dòng)量沿著氣動(dòng)外形表面當(dāng)?shù)鼐W(wǎng)格點(diǎn)的外法線方向施加,因此第i個(gè)POD基中第j個(gè)網(wǎng)格點(diǎn)的位置變化量(下文統(tǒng)稱為r-變量)為

其中:n為該網(wǎng)格點(diǎn)外法線方向.下文的設(shè)計(jì)知識(shí)主要是通過(guò)對(duì)r-變量進(jìn)行分析獲取.

3 氣動(dòng)外形優(yōu)化設(shè)計(jì)

3.1 同步并行粒子群方法

PSPSO方法是縮短基本PSO方法[9]優(yōu)化周期的有效手段.圖1為算法流程圖,其中n為粒子數(shù)目,m為優(yōu)化步數(shù),i為粒子編號(hào),k為當(dāng)前迭代步,f和c分別為目標(biāo)氣動(dòng)函數(shù)和約束函數(shù),pi為第i個(gè)粒子進(jìn)化到第k代時(shí)的最優(yōu)粒子,為第k代中所有粒子中的最優(yōu)粒子.粒子位置x和速度v的初始化方法與基本PSO相同,計(jì)算慣性權(quán)重因子采用基于差異系數(shù)(coefficient of variation,COV)的衰減方法[17].

圖1 PSPSO優(yōu)化流程圖Fig.1.Flowchart of PSPSO-based optimization.

3.2 參數(shù)化和網(wǎng)格變形方法

本文參數(shù)化方法是將葉片分解為多個(gè)特征截面,采用Hicks-Henne型函數(shù)方法[18]對(duì)特征截面進(jìn)行參數(shù)化,最后再通過(guò)插值得到整個(gè)三維葉片外形.本文參數(shù)化是在初始外形上疊加幾何擾動(dòng)量以產(chǎn)生新的外形,因此無(wú)需擬合初始外形,避免了在擬合初始外形時(shí)引入誤差.網(wǎng)格變形采用代數(shù)插值方法,將葉片表面的變化映射到空間網(wǎng)格,得到新外形的網(wǎng)格.為了提高計(jì)算效率,網(wǎng)格變形計(jì)算只對(duì)葉片相鄰的網(wǎng)格塊進(jìn)行.具體的參數(shù)化和網(wǎng)格變形方法見(jiàn)文獻(xiàn)[16].

4 Rotor 37數(shù)據(jù)挖掘算例

4.1 CFD數(shù)值驗(yàn)證及優(yōu)化問(wèn)題描述

流場(chǎng)計(jì)算采用團(tuán)隊(duì)開(kāi)發(fā)的大規(guī)模并行計(jì)算程序ASPAC,采用多塊結(jié)構(gòu)化網(wǎng)格.該程序求解雷諾平均Navier-Stokes方程,采用固連在轉(zhuǎn)子上的旋轉(zhuǎn)直角坐標(biāo)系來(lái)描述流動(dòng).湍流模型選用一方程Spalart-Allmaras模型[19].無(wú)黏通量采用Roe格式離散,黏性項(xiàng)采用二階中心差分格式離散.本文選擇跨音速壓氣機(jī)轉(zhuǎn)子葉片NASA Rotor 37進(jìn)行研究,該模型試驗(yàn)數(shù)據(jù)比較詳細(xì),是葉輪機(jī)CFD領(lǐng)域常用的驗(yàn)證算例之一,試驗(yàn)測(cè)量站位、空間分布及結(jié)果見(jiàn)文獻(xiàn)[20],設(shè)計(jì)參數(shù)見(jiàn)表1.

表1 Rotor 37設(shè)計(jì)參數(shù)Table 1.Design parameters of Rotor 37.

在出現(xiàn)旋轉(zhuǎn)失速之前,轉(zhuǎn)子葉片流場(chǎng)滿足周期性邊界條件,因此計(jì)算域僅選取一個(gè)葉片通道.圖2為O4H拓?fù)涞挠?jì)算網(wǎng)格.葉高方向分布56個(gè)網(wǎng)格單元,葉頂間隙內(nèi)徑向網(wǎng)格單元為16.流向和周向網(wǎng)格單元分別為152和60,總網(wǎng)格單元約為90萬(wàn).

圖2 Rotor 37計(jì)算網(wǎng)格 (a)S1流面;(b)S2流面Fig.2.Grid of Rotor 37:(a)In S1 streamsurface;(b)in S2 streamsurface.

圖3 (網(wǎng)刊彩色)Rotor 37特征工作線 (a)總壓比;(b)絕熱效率Fig.3. (color online)Performance characteristics:(a)Total pressure ratio;(b)adiabatic efficiency.

研究中還采用商業(yè)軟件NUMECA來(lái)進(jìn)行流場(chǎng)數(shù)值模擬.ASPAC堵塞流量為20.95 kg/s,NUMECA堵塞流量為20.93 kg/s,均在試驗(yàn)結(jié)果的誤差范圍內(nèi).圖3對(duì)比了CFD和試驗(yàn)測(cè)量的總壓、效率特征工作線.由圖3可知:基于CFD計(jì)算的總壓比均與試驗(yàn)較為接近,而CFD效率均與試驗(yàn)存在較明顯的偏差.截至目前,很少有數(shù)值模擬的效率與試驗(yàn)結(jié)果符合較好.Denton在[21]對(duì)Rotor 37的“盲測(cè)”結(jié)果進(jìn)行分析時(shí)指出:試驗(yàn)中端壁附近流動(dòng)損失無(wú)法測(cè)量可能是造成這種偏差的主要原因.通過(guò)與NUMECA計(jì)算結(jié)果的對(duì)比分析可知:ASPAC特征工作線稍優(yōu)于NUMECA,因而基于ASPAC的計(jì)算結(jié)果是合理的.另一方面,對(duì)于氣動(dòng)外形優(yōu)化設(shè)計(jì),更為重要的是能夠正確反映出氣動(dòng)性能隨設(shè)計(jì)變量的變化趨勢(shì).

4.2 基于粒子群的優(yōu)化設(shè)計(jì)

本文在Rotor 37的設(shè)計(jì)轉(zhuǎn)速近效率峰處進(jìn)行外形優(yōu)化設(shè)計(jì),同時(shí)約束近效率峰的總壓比和流量,優(yōu)化問(wèn)題可以由下式描述:

其中:η為等熵效率;m和m0分別為設(shè)計(jì)流量和目標(biāo)流量,目標(biāo)流量為原始外形的流量;π和π0分別為設(shè)計(jì)總壓比和目標(biāo)總壓比.

選擇位于100%,75%,50%,25%,0%葉高的五個(gè)特征截面進(jìn)行參數(shù)化,選擇中間三個(gè)截面進(jìn)行設(shè)計(jì).一般情況下,吸力面對(duì)絕熱效率的影響更大,只在吸力面布置10個(gè)型函數(shù),共30個(gè)設(shè)計(jì)變量.設(shè)計(jì)變量的變化范圍定義為(?0.1,0.1).粒子群優(yōu)化方法采用80個(gè)粒子,每個(gè)粒子使用8個(gè)CPU計(jì)算,一共使用640個(gè)CPU,共優(yōu)化25代.

表2給出了優(yōu)化前后的氣動(dòng)性能參數(shù),絕熱效率由0.8588提高到0.8692,提高了1.21%,總壓比降低了0.46%,流量提高了0.43%,均在約束條件之內(nèi).

表2 優(yōu)化前后葉片的氣動(dòng)性能Table 2.Aerodynamic performance of reference and optimized blades.

圖4給出了葉片吸力面上的壓力分布對(duì)比,50%葉高以上激波明顯向葉片尾緣移動(dòng),從而影響激波誘導(dǎo)的分離.圖5給出了吸力面附近的流線分布,圖中藍(lán)色流線為優(yōu)化后的流線,由圖可知:優(yōu)化后的激波位置明顯后移,相應(yīng)的分離區(qū)域明顯減小.

圖5 (網(wǎng)刊彩色)葉片吸力面流線(相對(duì)速度):原始(紅色);優(yōu)化(藍(lán)色)Fig.5.(color online)Streamlines(relative velocity)on the suction surface of reference(red)and optimized(blue)blades.

圖6給出了25%,50%,75%葉高的氣動(dòng)外形和壓力分布.25%葉高處,從弦長(zhǎng)中部開(kāi)始葉片厚度開(kāi)始減小,在尾緣附近先增大后減小,整體上最大厚度后移,能延遲激波的產(chǎn)生.由壓力分布也可以看出:優(yōu)化后激波明顯后移.50%葉高處,在尾緣附近葉片厚度明顯增加,最大厚度后移,激波也隨之后移,且激波后移更多.75%葉高處,從弦長(zhǎng)中部開(kāi)始葉片厚度開(kāi)始增大,在尾緣附近先減小后增大,但變化非常小.相對(duì)于25%葉高,75%葉高的優(yōu)化葉型雖然沒(méi)有明顯的最大厚度后移,但仍能延遲激波的產(chǎn)生.

由上述分析可知,Rotor 37絕熱效率的提高主要源自激波位置的后移,以及吸力面流動(dòng)分離區(qū)域的減小,導(dǎo)致這種變化的主要原因是葉型吸力面外形變化.

圖6 (網(wǎng)刊彩色)不同葉高的氣動(dòng)外形和壓力分布 (a)25%葉高;(b)50%葉高;(c)75%葉高Fig.6.(color online)Aerodynamic shape and pressure distributions at(a)25%span,(b)50%span,(c)75%span.

4.3 設(shè)計(jì)知識(shí)挖掘

優(yōu)化共進(jìn)行了25步,加上初始化粒子,共有2080個(gè)粒子,其中氣動(dòng)外形奇異和氣動(dòng)特性不滿足優(yōu)化要求的粒子在進(jìn)行篩選時(shí)需要剔除.本文旨在獲取與優(yōu)化設(shè)計(jì)目標(biāo)相關(guān)的設(shè)計(jì)知識(shí),篩選數(shù)據(jù)的規(guī)則為:滿足優(yōu)化約束且絕熱效率提高.篩選出有效粒子1734個(gè),并按照絕熱效率由高至低進(jìn)行排列.圖7給出了所有粒子的分類情況,其中:＜0%表示不滿足規(guī)則.

對(duì)滿足規(guī)則的數(shù)據(jù)按照(4)式進(jìn)行預(yù)處理獲得快照,并由SVD確定1734個(gè)POD基.POD基及特征值的特點(diǎn)是:特征值大的POD基包含最明顯的特征,本文采用基的能量比重對(duì)有效基進(jìn)行篩選,由(6)和(7)式計(jì)算的前10個(gè)POD基的能量比重如圖8所示.由圖可知,前三個(gè)基分別占了94.62%,3.65%和1.09%的總能量,99.36%的能量都集中在前三個(gè)基上面,因而下文主要對(duì)前三個(gè)基進(jìn)行分析.

圖7 (網(wǎng)刊彩色)所有粒子的分類情況Fig.7.(color online)Classi fication of proceeding data.

圖8 前10個(gè)基的能量比重Fig.8.Energy ratio of the first ten POD modes.

POD基的分析對(duì)象為基的r-變量,r-變量的正負(fù)并不直接反映葉片外形厚度的增大或減小,還要考慮基系數(shù)的影響.本文將r-變量的正負(fù)描述為正向變化和負(fù)向變化.由上文所述,從基能夠獲得的設(shè)計(jì)知識(shí)為氣動(dòng)外形變化量的相對(duì)大小和所處區(qū)域.為了更好地描述外形變化量的位置,本文將葉片吸力面劃分為9個(gè)區(qū)域,劃分方法為:沿葉高選擇25%,50%,75%三個(gè)截面,以1,2,3命名;沿弦向?qū)⑦@三個(gè)截面劃分為三個(gè)大致相等的區(qū)域,以A,B,C命名.圖9給出了9個(gè)區(qū)域的r-變量分布.由圖9可知:外形變形量大的區(qū)域位于B1,C1和B2,C2區(qū)域;變形復(fù)雜的區(qū)域位于B3,C3區(qū)域;A區(qū)域整體變形都很小.需要注意的是:r-變量的比較分析都是基于相同的原始外形幾何尺寸.根據(jù)上述設(shè)計(jì)知識(shí)可知:對(duì)于此類優(yōu)化設(shè)計(jì)問(wèn)題,設(shè)計(jì)變量應(yīng)重點(diǎn)分布在B,C區(qū)域,設(shè)計(jì)變量的設(shè)計(jì)空間在截面1,2處應(yīng)較大,在截面3處應(yīng)較小.

接下來(lái)分析基系數(shù)的影響.由于前三個(gè)POD基所占能量已經(jīng)超過(guò)99%,本文用前三個(gè)基來(lái)描述目標(biāo)函數(shù)排序靠前的100個(gè)快照.圖10給出了快照的前三個(gè)基系數(shù).由圖可知:三個(gè)基系數(shù)的絕對(duì)值依次減小,與基所占能量的大小變化趨勢(shì)一致;此外,基系數(shù)都在某一固定值附近波動(dòng),該固定值應(yīng)該能夠反映基系數(shù)中隱含的主要設(shè)計(jì)規(guī)律.

將基系數(shù)分解為固定值和擾動(dòng)的迭加形式:

其中αm表示前100個(gè)快照的第m個(gè)基的系數(shù)組成的向量,為該向量的固定值,為該向量擾動(dòng)值.采用線性回歸的方法確定上述公式中的固定值,可得:

其中n表示第n個(gè)快照.

圖10中不同顏色的虛線對(duì)應(yīng)各個(gè)基系數(shù)的線性回歸結(jié)果.雖然每個(gè)快照的固定值略有不同,但是差別非常小,可以忽略,因而

此外,三組基系數(shù)的擾動(dòng)量最大值分別為:

相對(duì)于固定值,前兩個(gè)基系數(shù)擾動(dòng)量的最大值較小,可忽略其對(duì)基權(quán)重的影響.第三個(gè)基系數(shù)擾動(dòng)量與固定值量級(jí)相同,但是考慮第三個(gè)基的權(quán)重僅占總能量的1.09%,其對(duì)外形的影響很小,基本可以忽略第三個(gè)基系數(shù)的擾動(dòng)量的影響.因而,可以僅由POD基及基系數(shù)的固定值確定“擬合優(yōu)化外形”,如下式所示:

其中φcomb表示POD基組合,也即擬合優(yōu)化外形.

圖9 (網(wǎng)刊彩色)r-變量分布云圖 (a)第一個(gè)基;(b)第二個(gè)基;(c)第三個(gè)基Fig.9.(color online)Contours of r-variation on:(a)The 1st mode;(b)the 2nd mode;(c)the 3rd mode.

圖10 (網(wǎng)刊彩色)前三個(gè)基的系數(shù)分布Fig.10.(color online)Coefficient distributions of the if rst three basis modes.

圖11給出了擬合優(yōu)化外形及優(yōu)化外形的r-變量分布.由圖可知:截面1,2,3上兩者的變化趨勢(shì)基本相同;截面1上,C1區(qū)域中部氣動(dòng)外形變化的最大值略有差別;截面2上,A2區(qū)域優(yōu)化外形減薄區(qū)域較擬合優(yōu)化外形的小,氣動(dòng)外形變化的最大值雖有差別,但都接近零,C2區(qū)域氣動(dòng)外形增厚的最大值略有差別;截面3上,B3區(qū)域氣動(dòng)外形增厚的最大值略有差別,C3區(qū)域氣動(dòng)外形增厚、減薄的最大值也略有差別.

由前三個(gè)POD基和基系數(shù)的固定值線性擬合得到的r-變量的變化規(guī)律與真實(shí)優(yōu)化外形的r-變量變化規(guī)律基本一致,對(duì)優(yōu)化結(jié)果和數(shù)據(jù)挖掘結(jié)果起到了相互驗(yàn)證的作用.

5 結(jié) 論

1)數(shù)據(jù)挖掘得到的設(shè)計(jì)知識(shí)與優(yōu)化設(shè)計(jì)結(jié)果基本一致,證實(shí)了本文采用的數(shù)據(jù)挖掘方法和優(yōu)化設(shè)計(jì)結(jié)果均是可靠的,驗(yàn)證了本文所采用的優(yōu)化設(shè)計(jì)方法和CFD程序;

2)基于POD的數(shù)據(jù)挖掘方法所獲取的設(shè)計(jì)知識(shí)能夠直觀反映氣動(dòng)外形的變化規(guī)律,具體到本文的葉片算例,明確地反映了設(shè)計(jì)變量的布置位置和變化幅度,同類型葉片均可按此規(guī)律進(jìn)行氣動(dòng)外形的優(yōu)化設(shè)計(jì),具有較大的指導(dǎo)意義;

3)數(shù)據(jù)挖掘的目標(biāo)還可以進(jìn)一步擴(kuò)展,不局限于優(yōu)化設(shè)計(jì)目標(biāo),比如:可以從效率優(yōu)化過(guò)程數(shù)據(jù)中提取壓比增大的過(guò)程數(shù)據(jù),對(duì)壓比增大的氣動(dòng)外形進(jìn)行數(shù)據(jù)挖掘.整體上,提高數(shù)據(jù)挖掘目標(biāo)的多樣性,能夠進(jìn)一步發(fā)掘優(yōu)化過(guò)程數(shù)據(jù)的隱含價(jià)值.

[1]Jeong S,Shimoyama K 2011Proc.Inst.Mech.Eng.Part G:J.Aerosp.Eng.225 469

[2]Jeong S,Chiba K,Obayashi S 2005J.Aeros.Comp.Inf.Com.2 452

[3]Chiba K,Obayashi S 2008J.Spacecraft Rockets45 975

[4]Oyama A,Nonomura T,Fujii K 2010J.Aircraft47 1756

[5]Oyama A,Verburg P,Nonomura T,Harry W M,Fujii K 2010AIAA Paper2010 1140

[6]Holmes P,Lumley J L,Berkooz G 1997Q.J.Roy.Meteor.Soc.123 2500

[7]Guo Z D,Song L M,Li J,Li G J,Feng Z P 2015J.Propul.Technol.36 207(in Chinese)[郭振東,宋立明,李軍,李國(guó)君,豐鎮(zhèn)平2015推進(jìn)技術(shù)36 207]

[8]Wang W,Mo R,Zhang Y 2013Comput.Eng.Appl.49 11(in Chinese)[汪偉,莫蓉,張巖2013計(jì)算機(jī)工程與應(yīng)用49 11]

[9]Sirvoich L,Kirby M 1987Quart.Appl.Math.45 561

[10]Duan Y H,Cai J S,Li Y Z 2012AIAA J.50 968

[11]LeGresley P,Alonso J 2000AIAA Paper2000 2545

[12]Toal D J J,Bresslo ffN W,Keane A J,Holden C M E 2010AIAA J.48 916

[13]Ghoman S,Wang Z,Chen P,Kapania K 2012AIAA Paper2012 1808

[14]Luo J,Duan Y,Tang X,Liu F 2015ASME Paper2015 42876

[15]Luo J Q,Duan Y H,Xia Z H 2016Acta Phys.Sin.65 124702(in Chinese)[羅佳奇,段焰輝,夏振華2016物理學(xué)報(bào)65 124702]

[16]Duan Y H,Wu W H,Fan Z L,Chen T 2016ASME Paper2016 56861

[17]Venter G,Sobieszczanski-Sobieski J 2003AIAA J.41 1583

[18]Hicks R M,Henne P A 1987J.Aircraft15 407

[19]Spalart P R A,Allmaras S 1992AIAA Paper1992 0439

[20]Reid L,Moore R D 1978NASA TP1978 1337

[21]Denton J D 1998J.Therm.Sci.6 1

PACS:02.60.Pn,47.85.Gj,47.11.—j,29.85.CaDOI:10.7498/aps.66.220203

*Projects supported by The National Nature Science Foundation of China(Grant Nos.51676003,51206003).

?Corresponding author.E-mail:duanyanhui@foxmail.com

?Corresponding author.E-mail:jiaqil@pku.edu.cn

Proper orthogonal decomposition-based data mining of aerodynamic shape for design optimization?

Duan Yan-Hui1)?Wu Wen-Hua1)Fan Zhao-Lin1)Luo Jia-Qi2)?

1)(China Aerodynamic Research and Development Center,Computational Aerodynamics Research Institute,

Mianyang 621000,China)2)(College of Engineering,Peking University,Beijing 100871,China)

2 July 2017;revised manuscript

18 July 2017)

Global optimization methods are becoming more and more important in aerodynamic shape optimization.A large number of proceeding data will be generated during design optimization,from which the implicit but valuable design knowledge can be extracted.The design knowledge can then be used to help the designers to acquire the e ff ects of geometric variations on the aerodynamic performance changes.In this paper,we strive to extract the implicit design knowledge from proceeding data by a data mining method based on proper orthogonal decomposition(POD),by which the design knowledge more enriched and more visualized than those obtained from other data mining methods can be obtained.Proceeding data for data mining are ingathered from aerodynamic shape optimization of a transonic compressor rotor blade,NASA Rotor 37.The design optimization attempts to maximize the adiabatic efficiency of Rotor 37 under the operation condition near peak efficiency with the constrains of mass flow rate and total pressure ratio.The parallel synchronous particle swarm optimization method is employed to search for the optimization in the design space.The particles with improved adiabatic efficiency,while within the optimization constrain tolerances are picked up from the design optimization,which are then used for data mining.The geometric coordinates of the aerodynamic shape with respect to the ingathered particles are regarded as the snapshots.Then the POD modes of the aerodynamic shape can be obtained by singular value decomposition on the snapshots.The results show that the universal rules of geometry variations for the optimization maximizing the adiabatic efficiency of Rotor 37 can be directly visualized by the design knowledge extracted from the proceeding data by POD-based data mining technique.Furthermore,the optimization results are also veri fied by the design knowledge extracted by data mining.

data mining,proper orthogonal decomposition,aerodynamic optimization design,transonic

10.7498/aps.66.220203

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):51676003,51206003)資助的課題.

?通信作者.E-mail:duanyanhui@foxmail.com

?通信作者.E-mail:jiaqil@pku.edu.cn