張威, 王強(qiáng), 路嘉晨, 閻超,*

(1. 北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100191; 2. 中國航天空氣動力技術(shù)研究院, 北京 100074)

在飛行器的整個壽命周期中存在著各種各樣的不確定性,而考慮不確定因素的氣動優(yōu)化研究在飛行器氣動外形設(shè)計中扮演著至關(guān)重要的角色。 考慮不確定性因素的優(yōu)化設(shè)計主要包括穩(wěn)健優(yōu)化設(shè)計[1]及可靠性優(yōu)化設(shè)計[2]。 穩(wěn)健設(shè)計的優(yōu)化目標(biāo)是為了在性能提升的同時降低目標(biāo)函數(shù)對不確定因素的敏感度,而可靠性設(shè)計則是為了降低故障概率。 本文的關(guān)注點(diǎn)主要在于穩(wěn)健優(yōu)化設(shè)計研究。 現(xiàn)有的穩(wěn)健優(yōu)化設(shè)計主要著重于針對環(huán)境擾動下的優(yōu)化設(shè)計研究,而對于幾何擾動的關(guān)注則相對較少。

對于幾何不確定性的考量,主要有2 方面:①局部位置的凸起凹陷,主要是由于工作狀態(tài)下的磨損、侵蝕等導(dǎo)致的。 由于其大小、位置、形狀、分布形式等較為依賴對大量實(shí)驗(yàn)數(shù)據(jù)的分析,而實(shí)驗(yàn)數(shù)據(jù)難以獲取使得該類幾何不確定性量化也較難實(shí)現(xiàn)。 ②加工過程中存在不確定性。 現(xiàn)有對該類幾何不確定性的研究中,由于認(rèn)知不足,或是根據(jù)經(jīng)驗(yàn)直接在如翼型前緣等關(guān)鍵位置添加滿足正態(tài)分布的幾何擾動[3],然后進(jìn)行不確定性量化;或是采用區(qū)間法[4]直接進(jìn)行幾何不確定性量化。 這些方法的使用,無論是對正態(tài)擾動添加位置的選擇與確定,還是區(qū)間法的直接應(yīng)用,經(jīng)驗(yàn)性都較強(qiáng),忽視了對于全局的考量,及潛在幾何模態(tài)的影響。

本文以RAE2822 翼型為例,采用主成分分析(principal component analysis, PCA)方法[5-6],對整個參數(shù)化過程進(jìn)行分析,增強(qiáng)對全局范圍內(nèi)主要幾何變形模態(tài)的認(rèn)知,了解對應(yīng)參數(shù)的分布規(guī)律。 應(yīng)用敏感度分析方法,針對不同目標(biāo)函數(shù)(以升力系數(shù)、阻力系數(shù)、翼型厚度等為例),找出主要變形控制模態(tài),作為幾何擾動來源。 同時,進(jìn)行了考慮幾何不確定性的穩(wěn)健優(yōu)化設(shè)計研究。

1 PCA-HicksHenne 方法

1.1 Hicks-Henne 參數(shù)化方法

Hicks-Henne 方法[7]是常用的參數(shù)化方法,本文以此對RAE2822 翼型進(jìn)行參數(shù)化,具體的幾何外形表達(dá)式如下:

式中:y為變形后翼型的上(下)表面函數(shù);yorigin為初始翼型的上(下)表面函數(shù);fk(x)為擾動型函數(shù),ck為對應(yīng)參數(shù)。

原始的Hicks-Henne 方法中的擾動型函數(shù)定義如下:

研究表明[8-9],該方法存在后緣擾動不足的問題。 對此,本文將翼型后緣擾動函數(shù)做了式(3)的修改。 其導(dǎo)函數(shù)(式(4))在后緣處取值不為0,可以提供有效的后緣擾動。

1.2 PCA 方法

PCA 是一種常用的數(shù)據(jù)降維方法。 通過對高維特征數(shù)據(jù)進(jìn)行降維,僅保留重要特征,從而實(shí)現(xiàn)數(shù)據(jù)處理效率的提升。

基于最近重構(gòu)性或最大可分性[10],可將PCA的優(yōu)化目標(biāo)寫為

式中:X為翼型樣本幾何向量組成的矩陣;W為協(xié)方差矩陣XXT的特征向量矩陣;I為單位矩陣。

采用拉格朗日乘子法可得

式中:wi為矩陣W第i列的向量,也是XXT的特征向量;λi為對應(yīng)的特征值。

因此,僅需對協(xié)方差矩陣XXT進(jìn)行特征值分解,得到特征值如式(7)所示。 一般對所得到的特征值按照由大到小順序進(jìn)行排序,特征值越大,表明對應(yīng)的“成分”的貢獻(xiàn)越大。 隨后再選取前d′個特征值對應(yīng)的特征向量W*= (w1,w2,…,wd′),即為主成分分析的解。 特征值由大到小如下:

1.3 PCA-HicksHenne 參數(shù)化方法

采用PCA 方法,對經(jīng)由Hicks-Henne 方法得到的大量翼型進(jìn)行分析,提取目標(biāo)翼型變形的主要模態(tài),并以此為基礎(chǔ)對翼型進(jìn)行變形。 該方法可以提高對翼型變形內(nèi)在幾何模態(tài)的認(rèn)知,進(jìn)一步降低參數(shù)化所需參數(shù)維度,提高效率,這就是本文所采用的PCA-HicksHenne 方法。

以RAE2822 翼型為例,設(shè)定Hicks-Henne 方法中的型函數(shù)個數(shù)為21(上下翼面合計變量維度為42)。 較多的型函數(shù)個數(shù)可以賦予參數(shù)化方法更好的設(shè)計空間表達(dá)能力。 在給定設(shè)計空間內(nèi)隨機(jī)選取100 000 個樣本翼型數(shù)據(jù),如圖1 所示,并對其進(jìn)行PCA 分析。 給出前20 個主要變形模態(tài)所對應(yīng)的能量分布,如圖2 和圖3 所示。 可以看出,前8 個變形模態(tài)表達(dá)了98% 以上的變形信息,這表明僅采用前8 個模態(tài)就能夠很好地表征原設(shè)計空間中的翼型變形。

圖1 樣本翼型及RAE2822 翼型示意圖Fig.1 Schematic of sample airfoils and RAE2822 airfoil

圖2 主要幾何變形模態(tài)對應(yīng)模態(tài)能量分布Fig.2 Energy amplitude of main geometric transformation modes

圖3 主要幾何變形模態(tài)對應(yīng)累積能量分布Fig.3 Accumulated energy of main geometric transformation modes

圖4 給出經(jīng)PCA 分析得到的前8 個變形模態(tài)。 模態(tài)1 和模態(tài)2 分別為翼型的厚度變形模態(tài)及彎度變形模態(tài),模態(tài)3 和模態(tài)4 分別為上下翼面的最大厚度位置的軸向位移模態(tài),模態(tài)5 ～模態(tài)8 則分別為翼面的擠壓模態(tài)。

圖4 主要幾何變形模態(tài)示意圖Fig.4 Schematic of main geometric transformation modes

本文以這些模態(tài)為設(shè)計變形模態(tài),通過線性疊加到基準(zhǔn)翼型RAE2822 上,從而實(shí)現(xiàn)翼型的變形。 同時,觀察本節(jié)隨機(jī)選取的100 000 個翼型對應(yīng)在這8 個模態(tài)構(gòu)建的參數(shù)空間中的參數(shù)分布情況,如圖5 所示。 觀察可知,所有模態(tài)的參數(shù)分布與N(0,0. 012)的正態(tài)分布幾乎完全一致,因此,可以認(rèn)為第3 節(jié)中所用到的不確定模態(tài)也滿足正態(tài)分布。 同時將各模態(tài)參數(shù)的設(shè)計范圍定義為[ -0.03,0.03],以滿足99.7% 以上的參數(shù)變化范圍。

圖5 主要幾何變形模態(tài)對應(yīng)參數(shù)分布情況Fig.5 Parameters distribution of main geometric transformation modes

2 重要參數(shù)識別

本文的主要目的是為了發(fā)展基于PCA-Hicks-Henne 參數(shù)化方法的幾何不確定性穩(wěn)健優(yōu)化設(shè)計方法。 因此,為了選定特定參數(shù)模態(tài)作為幾何擾動模態(tài),需要對各模態(tài)針對目標(biāo)函數(shù)(如升力系數(shù)、阻力系數(shù)、翼型厚度等)的參數(shù)敏感度進(jìn)行分析。

Sobol 指數(shù)方法是一種常用的全局參數(shù)敏感度分析方法[11-12]。 通過Sobol 指數(shù)方法,可以得到各參數(shù)對于目標(biāo)函數(shù)的影響,主要包括參數(shù)的主指標(biāo)、參數(shù)之間的交互指標(biāo)及參數(shù)的總指標(biāo)。參數(shù)的主指標(biāo)是指完全由該參數(shù)產(chǎn)生的影響;交互指標(biāo)是指參數(shù)之間相互作用,共同產(chǎn)生的影響,通常用來判斷參數(shù)之間是否存在關(guān)聯(lián);總指標(biāo)則是由參數(shù)主指標(biāo)與交互指標(biāo)共同組成,通常表示該參數(shù)產(chǎn)生的所有影響的總和。

本文以RAE2822 翼型的升力系數(shù)、阻力系數(shù)、翼型厚度為目標(biāo),研究第1 節(jié)中所得的各變形模態(tài)在設(shè)計空間內(nèi)的參數(shù)敏感度,結(jié)果分別如圖6 ～圖8 所示,計算狀態(tài)如表1 所示。 可以看出,對于升力而言,模態(tài)2 是主要影響模態(tài),這是由于跨聲速條件下翼型彎度是升力產(chǎn)生的主要來源;對于阻力而言,模態(tài)3 是主要影響模態(tài),這是由于在模態(tài)3 下,上翼面的最大厚度位置在軸向上發(fā)生變化,導(dǎo)致上翼面激波產(chǎn)生位移,激波強(qiáng)度、激波阻力發(fā)生變化,從而阻力發(fā)生變化;對于翼型厚度而言,模態(tài)1 即厚度變形模態(tài)是主要影響模態(tài)。

圖6 升力系數(shù)的參數(shù)敏感度分析Fig.6 Parametric sensitivity analysis of lift coefficient

圖7 阻力系數(shù)的參數(shù)敏感度分析Fig.7 Parametric sensitivity analysis of drag coefficient

圖8 翼型厚度的參數(shù)敏感度分析Fig.8 Parametric sensitivity analysis of airfoil thickness

表1 計算狀態(tài)Table 1 Computational condition

因此,模態(tài)1 ～模態(tài)3 是翼型性能的主要影響模態(tài)。 這些模態(tài)的微小變化,會對翼型性能產(chǎn)生較大的影響。 因此,本文選取這3 個模態(tài)作為主要幾何擾動模態(tài),進(jìn)行穩(wěn)健優(yōu)化設(shè)計。

3 幾何擾動定義

在進(jìn)行穩(wěn)健優(yōu)化設(shè)計之前,還需對幾何擾動模態(tài)的參數(shù)進(jìn)行定義。 根據(jù)第1 節(jié)研究結(jié)果,可以假定3 個擾動模態(tài)仍滿足正態(tài)分布,其均值為0。 為更好展現(xiàn)本文方法的性能,分別設(shè)定擾動模態(tài)參數(shù)的標(biāo)準(zhǔn)差為設(shè)計參數(shù)的0.2 ～0.4 倍,在基準(zhǔn)RAE2822 翼型上,探索這些不同標(biāo)準(zhǔn)差假設(shè)下幾何擾動的影響。

混沌多項(xiàng)式方法是一種廣泛使用的不確定度量化分析方法[13-14]。 本文選用二階混沌多項(xiàng)式方法來量化不同擾動量下的翼型升力系數(shù)、阻力系數(shù)、翼型厚度的不確定度,結(jié)果如表2 ～表4 所示。 可以看出,隨著幾何擾動量的增大,各個目標(biāo)函數(shù)的不確定性都在增強(qiáng)。 同時也注意到,0.4 倍擾動量下,99.7% 置信區(qū)間的翼型厚度變化范圍達(dá)到10 倍風(fēng)洞模型允許公差[15],遠(yuǎn)遠(yuǎn)偏離實(shí)際允許情況,不予考慮。

表2 不同擾動量下升力系數(shù)的不確定度對比Table 2 Comparisons of lift coefficient uncertainty under different perturbations

表3 不同擾動量下阻力系數(shù)的不確定度對比Table 3 Comparisons of drag coefficient uncertainty under different perturbations

表4 不同擾動量下翼型厚度的不確定度對比Table 4 Comparisons of airfoil thickness uncertainty under different perturbations

因此,在第4 節(jié)的穩(wěn)健優(yōu)化設(shè)計中,定義模態(tài)1 ～模態(tài)3 為幾何擾動模態(tài),符合正態(tài)分布,均值為0,標(biāo)準(zhǔn)差為0. 3 倍的設(shè)計參數(shù)標(biāo)準(zhǔn)差(即0.003)。 擾動下的翼型變化情況如圖9 所示。

圖9 幾何擾動下的翼型變化情況Fig.9 Airfoil transformation under geometric perturbation

4 穩(wěn)健優(yōu)化設(shè)計

4.1 優(yōu)化流程及計算說明

圖10 給出了幾何擾動下的穩(wěn)健優(yōu)化設(shè)計流程。 在無擾動的常規(guī)優(yōu)化流程中,先對初始外形進(jìn)行PCA-HicksHenne 參數(shù)化設(shè)計;其次對參數(shù)進(jìn)行試驗(yàn)設(shè)計,得到采樣點(diǎn);再次進(jìn)行配套網(wǎng)格變形及計算流體力學(xué)(computational fluid dynamic,CFD)解算,得到采樣點(diǎn)的響應(yīng);隨后在采樣點(diǎn)上構(gòu)建代理模型,并進(jìn)行優(yōu)化算法尋優(yōu);最終得到優(yōu)化外形。 而在考慮幾何擾動的優(yōu)化流程中,還需要在參數(shù)化完成之后進(jìn)行敏感度分析,尋找?guī)缀螖_動參數(shù),并根據(jù)CFD 結(jié)果進(jìn)行擾動下的不確定度量化,隨后再構(gòu)建代理模型、進(jìn)行算法尋優(yōu)。

圖10 幾何擾動下穩(wěn)健優(yōu)化設(shè)計流程Fig.10 Flow chart of robust optimization design under geometric perturbation

本文選取RAE2822 翼型作為基礎(chǔ)翼型,在自研CFD 計算平臺MI-CFD[16]上,構(gòu)建克里金模型[17],并利用多島遺傳算法[18]進(jìn)行尋優(yōu)。 計算狀態(tài)如表1 所示。 為展示穩(wěn)健優(yōu)化的性能,將其與確定優(yōu)化的結(jié)果進(jìn)行對比,具體見4.2 節(jié)和4.3 節(jié)。

4.2 案例1

本節(jié)案例是單目標(biāo)優(yōu)化案例。 確定優(yōu)化目標(biāo)及對應(yīng)的穩(wěn)健優(yōu)化目標(biāo)定義分別為

maxCl

式中:Cl和θ分別為翼型的升力系數(shù)和翼型厚度;μ和σ分別為對應(yīng)下標(biāo)變量在擾動下的均值和標(biāo)準(zhǔn)差。

優(yōu)化結(jié)果如表5 所示。 本文中無擾動翼型的升阻力性能與擾動下平均性能接近。 可以看出,確定優(yōu)化外形在升力均值方面提升很多,但同時也不可避免地伴隨著升力標(biāo)準(zhǔn)差的增長。 性能提升的同時,不確定度也在增大,這是不利的。 而對于穩(wěn)健優(yōu)化外形,因?yàn)樵趦?yōu)化過程中已經(jīng)考慮了對不確定性因素的約束,所以在其性能提升的同時,不確定度基本無變化。

表5 案例1 優(yōu)化結(jié)果對比Table 5 Comparisons of optimal results in Case 1

4.3 案例2

本節(jié)案例是多目標(biāo)優(yōu)化案例。 確定優(yōu)化目標(biāo)及穩(wěn)健優(yōu)化目標(biāo)定義分別為

式中:Cd為翼型的阻力系數(shù)。

圖11 和圖12 給出2 個優(yōu)化過程的Pareto前沿面,其中穩(wěn)健優(yōu)化結(jié)果僅給出平均升力與平均阻力的對應(yīng)關(guān)系。 在圖11 和圖12 中分別選出3 個設(shè)計點(diǎn),并將結(jié)果在表6 和表7 中列出。

圖11 案例2 確定優(yōu)化Pareto 前沿面Fig.11 Pareto frontier of deterministic optimization in Case 2

圖12 案例2 穩(wěn)健優(yōu)化Pareto 前沿面Fig.12 Pareto frontier of robust optimization in Case 2

表6 案例2 確定優(yōu)化結(jié)果Table 6 Solutions of deterministic optimization in Case 2

表7 案例2 穩(wěn)健優(yōu)化結(jié)果Table 7 Solutions of robust optimization in Case 2

如表6 和表7 所示,所有得到的優(yōu)化翼型均滿足翼型厚度約束。 其中,確定優(yōu)化得到的3 個優(yōu)化翼型在平均升阻力性能得到提升的同時,對阻力系數(shù)的不確定度恰好也有減弱的效果,而對升力系數(shù)的不確定度則恰恰相反。 3 個翼型的升力系數(shù)標(biāo)準(zhǔn)差均有不同幅度的增長,最高達(dá)到了近200%。 對于飛行器穩(wěn)健性的設(shè)計要求,這是不利的。 而穩(wěn)健優(yōu)化得到的3 個優(yōu)化翼型,則在目標(biāo)性能相斥的情況下,兼顧了升阻力的均值及標(biāo)準(zhǔn)差等各方面的性能要求。 在提升某些方面性能的同時,也能保證其他性能的不惡化。

由表6 和表7 確定優(yōu)化得到的翼型③(以下為確定優(yōu)化翼型)和穩(wěn)健優(yōu)化得到的翼型⑤(以下稱穩(wěn)健優(yōu)化翼型)的平均性能相近,對其作對比分析。

圖13 給出了RAE2822 翼型、確定優(yōu)化翼型及穩(wěn)健優(yōu)化翼型的形狀對比,并在圖14 ～圖16中給出了3 個翼型及其對應(yīng)擾動翼型的壓力系數(shù)Cp的分布情況。

圖13 優(yōu)化翼型對比Fig.13 Comparisons of optimal airfoils

圖14 RAE2822 翼型及擾動翼型的壓力分布Fig.14 Surface pressure distribution of RAE2822 airfoil and its corresponding airfoils under perturbation

圖15 確定優(yōu)化翼型及擾動翼型的壓力分布Fig.15 Surface pressure distribution of deterministic optimal airfoil and its corresponding airfoils under perturbation

如圖13 ～圖16 所示,在設(shè)計狀態(tài)下,2 個優(yōu)化外形在下翼面后緣附近都存在由于彎度增加的后緣加載,升力因此增加。 此外,穩(wěn)健外形在上翼面前緣附近曲率更大,產(chǎn)生了一個較小的吸力峰,因此其升力相對更大。 而在阻力方面,2 個優(yōu)化外形的上翼面相對RAE2822 翼型都更加平緩,確定優(yōu)化外形分解為2 個小激波,穩(wěn)健優(yōu)化外形激波強(qiáng)度減弱,二者阻力均降低。

圖16 穩(wěn)健優(yōu)化翼型及擾動翼型的壓力分布Fig.16 Surface pressure distribution of robust optimal airfoil and its corresponding airfoils under perturbation

而在幾何擾動下,穩(wěn)健優(yōu)化外形的激波強(qiáng)度和位置的變化與RAE2822 翼型相似,因此兩者的升力系數(shù)變化幾乎一致,標(biāo)準(zhǔn)差相當(dāng)。 另外穩(wěn)健外形的激波強(qiáng)度相對較小,因此其阻力系數(shù)的標(biāo)準(zhǔn)差也相對較小。 對于確定優(yōu)化外形,其激波位置在擾動下變化劇烈,由此導(dǎo)致升力變化較大。而其激波強(qiáng)度較弱,變化也不大,因此阻力變化不大。

5 結(jié) 論

本文構(gòu)建了基于主成分分析的PCA-HicksHenne參數(shù)化方法,并在此基礎(chǔ)上進(jìn)行了考慮幾何不確定性的穩(wěn)健優(yōu)化研究,結(jié)論如下:

1) 對原始Hicks-Henne 方法進(jìn)行了修改,使其在翼型后緣處也能提供有效擾動,提高了參數(shù)化方法的空間表達(dá)能力。

2) 構(gòu)建了基于主成分分析的PCA-HicksHenne參數(shù)化方法,揭示了翼型的主要幾何變形模態(tài),驗(yàn)證了對應(yīng)參數(shù)的分布規(guī)律,同時實(shí)現(xiàn)了參數(shù)化方法的降維,提高了優(yōu)化效率。

3) 利用參數(shù)敏感度分析方法,發(fā)現(xiàn)翼型厚度變形模態(tài)、彎度變形模態(tài)及上翼面最大厚度位置的軸向位移模態(tài)是翼型性能的主要影響模態(tài),必須予以重視。

4) 發(fā)展了基于PCA-HicksHenne 方法幾何不確定性穩(wěn)健優(yōu)化設(shè)計方法,拓展了對幾何不確定性穩(wěn)健優(yōu)化的研究,優(yōu)化得到的目標(biāo)翼型在性能提升的同時,穩(wěn)健性也得以保證。