(西南科技大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，四川 綿陽 621010)

0 引言

近幾十來世界航空航天技術(shù)取得了迅猛的發(fā)展，航空航天已成為大國間的主要競技場之一，也是國家綜合實力的體現(xiàn)，而飛行器的設(shè)計又是航空航天的重要研發(fā)領(lǐng)域。機翼的設(shè)計是整個飛行器設(shè)計的核心之一，而翼型又是機翼設(shè)計的核心。如何快速、低成本地獲得一個性能高效且穩(wěn)定的翼型成為了研究人員關(guān)注的重點之一。

翼型的傳統(tǒng)設(shè)計方法主要是通過風(fēng)洞對研究人員建立的不同機翼模型進行試驗，通過試驗數(shù)據(jù)分析來選擇最佳翼型，該種方法存在成本高、試驗周期長等缺點。此外傳統(tǒng)翼型設(shè)計中，設(shè)計人員在借鑒其過去設(shè)計經(jīng)驗時，可能會引起單點設(shè)計問題如對單一點馬赫數(shù)處的阻力進行了優(yōu)化設(shè)計，但引發(fā)了該點臨域馬赫數(shù)處的阻力出現(xiàn)劇烈波動的情況。在超聲速翼型的設(shè)計中單點設(shè)計帶來的性能波動會更加劇烈，在Hicks和Johnson的研究[1]中對單點設(shè)計問題進行了詳細介紹。

在之后的研究中，研究人員們將計算流體動力學(xué)(Computational Fluid Dynamics，CFD)等相關(guān)數(shù)值計算用于翼型穩(wěn)健設(shè)計。雖然通過CFD軟件進行氣動性能計算能降低傳統(tǒng)風(fēng)洞試驗的成本，但是其計算時間仍然非常漫長。而通過將替代模型與CFD計算相結(jié)合的方法來進行翼型穩(wěn)健設(shè)計，能夠在保證精度需求的情況下大大縮減計算時間，能在基準(zhǔn)翼型的基礎(chǔ)上獲得一個性能更優(yōu)的穩(wěn)健翼型。

1 翼型穩(wěn)健設(shè)計及其目標(biāo)模型

翼型的穩(wěn)健設(shè)計就是要對選定的基準(zhǔn)翼型進行優(yōu)化，以最終獲得一個對環(huán)境因素變化不敏感的穩(wěn)健翼型。因此在翼型設(shè)計的初始階段就要明確穩(wěn)健設(shè)計針對的噪聲因素，通過穩(wěn)健設(shè)計方法來獲取一組可控因素的最佳組合，從而實現(xiàn)穩(wěn)健翼型對噪聲因素不敏感，達到翼型性能高且穩(wěn)定的目的。

飛行器翼型在非設(shè)計狀態(tài)下常出現(xiàn)性能不穩(wěn)定現(xiàn)象。本文以一馬赫數(shù)區(qū)間為非設(shè)計狀態(tài)因素，進行RAE2822基準(zhǔn)翼型的穩(wěn)健設(shè)計，使穩(wěn)健翼型的阻力系數(shù)對馬赫數(shù)區(qū)間范圍內(nèi)的變化不敏感。達到不但要降低翼型的阻力還要實現(xiàn)翼型性能穩(wěn)定即不發(fā)生阻力系數(shù)劇烈波動的目的。該目的在本文中通過對穩(wěn)健翼型在馬赫數(shù)區(qū)間(Ma∈[MaminMamax])的阻力系數(shù)的均值和方差進行評價，本文穩(wěn)健設(shè)計參考了文獻2[2]中的目標(biāo)模型，本文目標(biāo)模型如下所示：

(1)

式中，μ為阻力系數(shù)的均值；σ2為阻力系數(shù)的方差；D為翼型的幾何外形參數(shù)；T為外形的幾何約束；Cl為升力系數(shù)約束條件；R為雷諾數(shù)約束條件；α為攻角；Ma為選定的馬赫數(shù)變化區(qū)間。

在穩(wěn)健設(shè)計目標(biāo)模型中翼型的升力系數(shù)和阻力系數(shù)都是關(guān)于翼型外性設(shè)計變量、攻角以及馬赫數(shù)的函數(shù)。穩(wěn)健設(shè)計的目標(biāo)是獲得一個在給定的馬赫數(shù)變化區(qū)間內(nèi)，較基準(zhǔn)翼型的阻力系數(shù)均值和方差更小的穩(wěn)健翼型。將升力系數(shù)和雷諾數(shù)作為約束條件，避免其對最終阻力系數(shù)結(jié)果產(chǎn)生影響。使用偏最小二乘法(Partial Least Squares，PLS)建立阻力系數(shù)關(guān)于翼型外形設(shè)計參數(shù)(D)和馬赫數(shù)(Ma)的替代模型來對阻力系數(shù)進行預(yù)測。PLS替代模型可解決使用CFD軟件計算阻力系數(shù)Cd時計算量過大、耗時長，甚至無法求解的問題,并能保證最終結(jié)果滿足一定的精度要求。使用PLS替代模型的預(yù)測結(jié)果來獲得目標(biāo)模型≠2+...2的近似值，目標(biāo)模型的近似值求解公式如下所示：

(2)

2 翼型穩(wěn)健設(shè)計流程及各步驟原理

本文中使用偏最小二乘法建立替代模型進行翼型穩(wěn)健的流程主要分為三大關(guān)鍵步驟：確定試驗設(shè)計方法、建立替代模型和進行遺傳優(yōu)化。穩(wěn)健設(shè)計詳細步驟與設(shè)計流程圖(圖1)如下所示：

(1)確定穩(wěn)健設(shè)計的優(yōu)化目標(biāo)模型。

(2)進行翼型參數(shù)化，確定外形設(shè)計因素的個數(shù)。

(3)確定試驗設(shè)計方法及試驗設(shè)計的因素數(shù)及其水平數(shù)，生成對應(yīng)樣本點數(shù)據(jù)。

(4)對試驗設(shè)計樣本點翼型外形生成對應(yīng)網(wǎng)格并進行雷諾數(shù)約束下的CFD計算，獲得樣本點對應(yīng)的阻力系數(shù)Cd值。

(5)根據(jù)CFD計算出的樣本Cd的結(jié)果，使用偏最小二乘法建立替代模型。

(6)使用偏最小二乘法建立的替代模型與遺傳算法結(jié)合進行尋優(yōu)，得到穩(wěn)健翼型外形設(shè)計參數(shù)水平值的最佳組合。

(7)計算穩(wěn)健翼型外形參數(shù)所對應(yīng)翼型在設(shè)定的馬赫數(shù)范圍內(nèi)的阻力系數(shù)的均值與方差，并與基準(zhǔn)翼型RAE2822進行比較，驗證穩(wěn)健翼型性能是否更佳。

圖1 PLS翼型穩(wěn)健設(shè)計流程圖

穩(wěn)健翼型設(shè)計流程中各步驟詳細介紹：

2.1 翼型外形設(shè)計參數(shù)化

在翼型外形參數(shù)化步驟中采用當(dāng)前常用的Hicks-Henne[3]型函數(shù)進行線性疊加的方法來獲取樣本翼型的外形，樣本翼型的外形是由基準(zhǔn)翼型和型函數(shù)及其相關(guān)系數(shù)疊加生成的，Hicks-Henne生成樣本翼型的公式如下所示：

(3)

式中，F(xiàn)bas為基準(zhǔn)翼型的外形數(shù)據(jù)，Pi為型函數(shù)的相關(guān)系數(shù)，即翼型外形參數(shù)的設(shè)計變量，n為選取的設(shè)計變量的個數(shù)，fi(x)為Hicks-Henne的型函數(shù)。Hick-Henne型函數(shù)的表達式如下：

(4)

2.2 試驗設(shè)計方法的選取

試驗設(shè)計[4](Design of Experiment)是用于對試驗進行合理安排，使樣本點更具代表性，減少試驗次數(shù)的方法。常用的試驗設(shè)計方法主要有：正交試驗設(shè)計方法、均勻設(shè)計法、超拉丁方抽樣法等。本文采用由我國方開泰教授和數(shù)學(xué)家王元提出的均勻設(shè)計法。均勻設(shè)計法與正交試驗法相比，只考慮了試驗點在試驗設(shè)計范圍內(nèi)均勻分布，無需考慮正交試驗中整齊可比的要求，其試驗點的均勻性更佳，將其用于翼型穩(wěn)健設(shè)計能大大減少樣本點的數(shù)量。在本文中，通過馬赫數(shù)和翼型外形參數(shù)的因素數(shù)及水平數(shù)來確定最合適的均勻設(shè)計表，來進行試驗樣本的安排。

2.3 樣本翼型的網(wǎng)格生成及其CFD計算

對樣本數(shù)據(jù)表中每一行的一個翼型樣本使用Gridgen軟件進行對應(yīng)的二維網(wǎng)格生成，網(wǎng)格結(jié)構(gòu)采用C型結(jié)構(gòu)。在進行翼型樣本網(wǎng)格生成前，需要先建立RAE2822基準(zhǔn)翼型的結(jié)構(gòu)化二維網(wǎng)格，為后續(xù)CFD結(jié)果正確性驗證做準(zhǔn)備。RAE2822基準(zhǔn)翼型的C型網(wǎng)格如圖2所示。

圖2 RAE2822基準(zhǔn)翼型網(wǎng)格

生成樣本點翼型的網(wǎng)格后，開始進行對應(yīng)網(wǎng)格的CFD[5]計算，使用氣動計算中常用的Fluent軟件完成。在樣本翼型CFD計算前，先使用Fluent計算RAE2822基準(zhǔn)翼型網(wǎng)格在雷諾數(shù)R6=6.5×106，Ma=0.73，α=3.19時的翼型表面壓力分布，通過和1979年關(guān)于RAE2822的試驗報告[6]中的翼型表面壓力分布結(jié)果進行對比，來驗證CFD計算的正確性。對比結(jié)果如圖3所示，結(jié)果表明氣動分析結(jié)果正確，符合精度要求。

圖3 RAE2822翼型表面壓力分布CFD計算和試驗報告結(jié)果對比

2.4 建立偏最小二乘替代模型

建立偏最小二乘替代模型的目的是為了降低使用Fluent軟件進行CFD計算時的工作量，減少求解時間。通過PLS進行回歸擬合，能最快得到阻力系數(shù)目標(biāo)模型值。偏最小二乘法又稱偏最小二乘回歸[7]，其是一種多對多的線形回歸建模方法。當(dāng)所求解目標(biāo)的自變量和因變量數(shù)量多，且彼此存在多重相關(guān)性，而樣本的數(shù)量卻較少時，用偏最小二乘法獲得的模型比傳統(tǒng)回歸模型，如主成分回歸分析模型要更優(yōu)，能提供更加深入和豐富的信息。

2.4.1 偏最小二乘法的建模原理：

假設(shè)有m個自變量{x1,x2…,xp}和n個因變量{y1,y2…,yp}。為了研究自變量x和因變量y之間的統(tǒng)計關(guān)系,設(shè)有k個樣本點,由此獲得了自變量與因變量的數(shù)據(jù)集X={x1,x2…,xp}和Y={y1,y2…,yp}。

偏最小二乘法首先分別在自變量集合X與因變量集合Y中分別提取出第一成分t1和u1(簡而言之，t1是x1,x2…,xp的一個線形組合,u1是y1,y2…,yq的一個線形組合)。在提取這兩個第一成分時,為了回歸分析的需要,有兩個需要格外關(guān)注的要求:

(1)t1和u1應(yīng)盡可能多地分別提取出自變量集和因變量集中的變異信息;

(2)t1和u1間的相關(guān)程度要達到最大。

上述兩個要求的目的是要使得t1和u1盡可能好的代表數(shù)據(jù)集合X和Y,同時自變量的第一成分t1對因變量的第一成分u1還要具有最強的解釋能力。

在第一成分t1和u1被提取后，使用偏最小二乘法分別進行X對t1的回歸以及Y對u1的回歸。如果回歸方程能達到滿意的精度，則算法中止；否則將利用X被t1解釋后的殘余信息以及Y被u1解釋后的殘余信息進行第二成分的提取，直到能達所需要的精度為止。若最終對自變量X共提取出r個成分t1,t2…,tr，偏最小二乘法將通過建立y1,y2…,yq與t1,t2…,tr的回歸,然后再表達成y1,y2…,yq關(guān)于原變量x1,x2…,xp的回歸方程。

2.4.2 偏最小二乘法的簡化算法步驟：

本文采用文獻[10]中偏最小二乘的簡化算法，能極大的減少計算時間，減小程序代碼的編寫難度。

為簡便計算，設(shè)m個自變量{x1,x2…,xp}和n個因變量{y1,y2…,yq}均為歸一化后的變量。歸一化方法采用min-max方法，如下所示：

(5)

式中，i為變量范圍內(nèi)的一個具體變量的下標(biāo)；j為變量類型下標(biāo)(自變量或因變量)xjmax；xjmin和分布為自變量或因變量的最大值和最小值。

完成歸一化后，將自變量集和因變量集的w次數(shù)據(jù)矩陣分別記為：

PLS簡記算法的簡要步驟如下：

(6)

最終的n個因變量的偏最小回歸方程為 ：

yj=aj1x1+…+ajmxm,j∈[1,m]

(7)

偏最小二乘法的詳細步驟、推導(dǎo)過程及其Matlab程序代碼可參考文獻[10]。偏最小二乘法還可建立二次、三次等高階多項式回歸模型。本文中為獲得更好的模型擬合效果，采用偏最小二乘法建立三次多項式回歸模型。

2.5 遺傳算法優(yōu)化

本文采用遺傳算法進行穩(wěn)健翼型的尋優(yōu)。采用馬赫數(shù)按等差數(shù)列生成，其余翼型設(shè)計變量在指定范圍內(nèi)隨機生成的辦法來進行遺傳算法優(yōu)化求解。

遺傳算法[11]是一種自適應(yīng)全局搜索算法。其參考達爾文的生物進化論，通過模擬自然界中生物種群的進化發(fā)展過程，來盡可能得到全局最優(yōu)解。

本文中對目標(biāo)模型進行近似表示，來進行遺傳算法求解：

(8)

式中，NMa為選取的馬赫數(shù)區(qū)間范圍內(nèi)的馬赫數(shù)的個數(shù)，D為翼型的外形參數(shù)變量。

在遺傳算法算法計算過程中，對多個變量使用二進制的形式進行編碼。遺傳算法中種群的數(shù)量決定了遺傳算法的多樣性，種群數(shù)目越多則多樣性越好，但計算量也會增加，降低運行效率。然后數(shù)目過少又可能出現(xiàn)局部最優(yōu)解，因此需要研究人員根據(jù)具體問題進行種群數(shù)目的確定，一般種群數(shù)目最小應(yīng)大于等于20。在遺傳算法中變異概率(Pm)對種群的影響應(yīng)遠遠小于交叉概率(Pc)，Pc和Pm的相關(guān)取值分析見文獻[12]，Pc一般取值范圍為0.4～0.99，Pm取值范圍為0.001～0.1。罰函數(shù)采用的是翼型最大相對厚度的幾何約束，即最大相對厚度的最大值和最小值，使用罰函數(shù)能保證最終翼型外形的合理性，以便獲得最優(yōu)解。將PLS替代模型和遺傳算法結(jié)合以獲得穩(wěn)健翼型的外形設(shè)計參數(shù)的一組最佳組合，再使用Hicks-Henne型函數(shù)進行處理，即可獲得穩(wěn)健翼型的外形數(shù)據(jù)。

3 試驗樣例及其結(jié)果分析

以對RAE2822基準(zhǔn)翼型進行穩(wěn)健設(shè)計來對偏最小二乘法建立的替代模型的效果進行驗證。選擇在馬赫數(shù)Ma∈[0.7,0.8]，雷諾數(shù)Re=6.5×106,翼型升力系數(shù)Cl=0.8，翼型最大相對厚度0.1≤d≤0.12的條件下，進行基準(zhǔn)翼型的穩(wěn)健設(shè)計。在Hicks-Henne翼型參數(shù)化時采用十個設(shè)計參數(shù)變量(d1～d10)進行翼型外形的表示(上下翼面各取五個參數(shù)變量并分別確定各設(shè)計變量的取值范圍，如d1∈[-0.006,0.006]等)。馬赫數(shù)是設(shè)計中除了翼型外形十個設(shè)計變量外的另一個設(shè)計變量。

完成樣本CFD計算后進行偏最小二乘替代模型的建立。在PLS替代模型建立前先求取PLS模型三次多項式各系數(shù)的值，并進行保存。PLS替代模型中取NMa=1 000，即按公差為0.000 1在Ma∈[0.7,0.8]，選擇種群規(guī)模M= 100,交叉概率Pc= 0.8,變異概率為pm= 0.1，進化代數(shù)步長Pd= 50，初始種群進化代數(shù)=Pd。當(dāng)某進化步長的整數(shù)倍時的最佳翼型的目標(biāo)模型值與臨近左右步長的最佳翼型的目標(biāo)模型值相等時，即目標(biāo)模型值不再發(fā)生變化時，停止種群進化迭代。將穩(wěn)健翼型的外形設(shè)計變量組合的值帶入Hicks-Henne中，獲得穩(wěn)健翼型的外形數(shù)據(jù)，并生成對應(yīng)的二維結(jié)構(gòu)化網(wǎng)格。按公差為 0.001在中生成100個馬赫數(shù)，使用Fluent計算穩(wěn)健翼型分別在這100個馬赫數(shù)以及升力系數(shù)、雷諾數(shù)約束下的Cd值，計算出阻力系數(shù)的均值與方差并與基準(zhǔn)翼型在同樣條件下的均值和方差進行對比，以驗證穩(wěn)健翼型的效果是否更優(yōu)。

圖4 遺傳算法進化次數(shù)及對應(yīng)目標(biāo)模型近似值

遺傳算法計算環(huán)境為個人PC，硬件環(huán)境：處理器 2.7 GHz Intel Core i7，內(nèi)存16GB 2133MHz LPDDR3；軟件環(huán)境 Python 3.7.1，PyCharm 2018。遺傳算法的進化次數(shù)的計算結(jié)果如圖4所示，可以發(fā)現(xiàn)在進化次數(shù)為250次時，第200次的目標(biāo)模型的近似值已與第150次和第250次的值相等，不再發(fā)生改變，已求得遺傳算法的最優(yōu)解，停止進化迭代。遺傳算法的最終進化次數(shù)250次，用時112分鐘30秒。遺傳算法的最終結(jié)果如表1，可以看到偏最小二乘法建立的替代模型與遺傳算法結(jié)合獲得的穩(wěn)健翼型在時，穩(wěn)健翼型阻力系數(shù)的均值與方差都明顯小于RAE2822基準(zhǔn)翼型的數(shù)值，說明使用偏最小二乘建模方法獲得的穩(wěn)健翼型的效果更佳，偏最小二乘法建立替代模型的方法確實可以用于翼型的穩(wěn)健設(shè)計之中。

表1 PLS穩(wěn)健翼型與基準(zhǔn)翼型阻力系數(shù)均值和方差對比

圖5為PLS穩(wěn)健翼型與RAE2822基準(zhǔn)翼型的外形對比，可明顯發(fā)現(xiàn)穩(wěn)健翼型的厚度有所減小，上表面更為平坦，因此其出現(xiàn)阻力發(fā)散的臨界馬赫數(shù)將提高，最大相對厚度的減小能帶來翼型波阻的減小，使性能得到提升。此外穩(wěn)健翼型下表面后緣出現(xiàn)了一個非常明顯的向里凹進去的反曲斷，表明穩(wěn)健翼型具有了明顯的超臨界翼型的特征，后緣反曲斷的出現(xiàn)使后緣的升力增加，彌補了由于上表面平坦而引發(fā)的升力不足問題。在圖6翼型性能比較圖中，可明顯發(fā)現(xiàn)PLS穩(wěn)健翼型在升力系數(shù)和雷諾數(shù)約束下且處于設(shè)定范圍時，其阻力系數(shù)明顯小于基準(zhǔn)RAE2822翼型的阻力系數(shù)值。從最終結(jié)果來看，PLS建立的替代模型達到了翼型穩(wěn)健設(shè)計的目標(biāo)。

圖5 PLS穩(wěn)健翼型與基準(zhǔn)翼型的外形對比

圖6 PLS穩(wěn)健翼型與基準(zhǔn)翼型的性能對比

4 結(jié)論

本文使用偏最小二乘法建立替代模型來解決傳統(tǒng)飛行器翼型設(shè)計中存在的單點設(shè)計和CFD計算量過大，風(fēng)洞試驗成本高的問題。以獲得某一區(qū)間變化的馬赫數(shù)(即外界環(huán)境不穩(wěn)定)下翼型阻力系數(shù)均值和方差最低的翼型穩(wěn)健設(shè)計為研究背景，通過偏最小二乘法建立替代模型來完成翼型的穩(wěn)健設(shè)計，實現(xiàn)在保證一定精度下，穩(wěn)健設(shè)計快捷和簡易的目標(biāo)。結(jié)果表明，采用PLS替代模型方法，能有效實現(xiàn)本文翼型穩(wěn)健設(shè)計計算成本低，計算快，精度能保證的目的。

在后續(xù)的研究中，還可將PLS替代模型方法用于機翼的穩(wěn)健設(shè)計以及風(fēng)力機葉片穩(wěn)健設(shè)計、發(fā)動機葉片穩(wěn)健設(shè)計中。飛行器翼型的穩(wěn)健設(shè)計結(jié)果可為后續(xù)機翼和整機的穩(wěn)健設(shè)計打下了良好的基礎(chǔ)。偏最小二乘替代模型在穩(wěn)健設(shè)計上具有廣闊的應(yīng)用前景，值得深入研究。