胡漢鐸，宋彥萍，俞建陽(yáng)，劉瑤，陳浮，高文秀

哈爾濱工業(yè)大學(xué) 能源科學(xué)與工程學(xué)院，哈爾濱 150001

隨著航空技術(shù)的不斷發(fā)展，對(duì)飛行器的飛行性能要求不斷提高［1］。機(jī)翼作為飛行器升力的主要提供部件，其性能水平直接決定了飛行器的氣動(dòng)性能和安全裕度［2］。然而，在實(shí)際加工、裝配和飛行過程中，由于制造誤差［3］、結(jié)冰［4-5］和腐蝕［6］等因素的存在，機(jī)翼的外形不可避免地與設(shè)計(jì)的理想形狀存在偏差，翼型幾何在一定的范圍內(nèi)隨機(jī)變化，即存在幾何不確定性。不確定性因素的存在會(huì)使機(jī)翼性能發(fā)生波動(dòng)，進(jìn)而導(dǎo)致飛行器性能惡化，甚至帶來(lái)安全隱患［7］，因此開展不確定性量化、獲得機(jī)翼性能參數(shù)在考慮不確定性因素下的隨機(jī)響應(yīng)特性，對(duì)于評(píng)估飛行器安全性能、指導(dǎo)機(jī)翼魯棒性設(shè)計(jì)具有重要研究意義［8］。

在不確定性量化方法方面，蒙特卡洛采樣（Monte Carlo Sampling，MCS）和混沌多項(xiàng)式展開（Polynomial Chaos Expansion，PCE）已經(jīng)被國(guó)內(nèi)外學(xué)者廣泛應(yīng)用于翼型的不確定性研究中。Loeven等［9］于2007年提出用概率配點(diǎn)法（Probabilistic Collocation Method，PCM）采樣求解混沌多項(xiàng)式系數(shù)，并將其用于NACA0012翼型的不確定性量化，只考慮來(lái)流馬赫數(shù)的變化，根據(jù)概率配點(diǎn)法，只需要5個(gè)配點(diǎn)即可量化翼型升力和阻力的概率密度分布；不久后，Loeven等［10］將概率配點(diǎn)法應(yīng)用于考慮3個(gè)幾何不確定性因素的NACA5412翼型不確定性量化，配點(diǎn)數(shù)增加到27個(gè)。高遠(yuǎn)等［11］使用3階非嵌入式混沌多項(xiàng)式展開，針對(duì)考慮來(lái)流馬赫數(shù)和攻角不確定性的RAE2822翼型開展不確定性量化，使用25個(gè)樣本可以在效率與精度之間取得平衡。于佳鑫等［12］針對(duì)S809翼型，考慮了12個(gè)變量的幾何不確定性，對(duì)比了蒙特卡洛采樣和稀疏網(wǎng)格配點(diǎn)法獲得收斂量化結(jié)果所需的樣本數(shù)和精度，對(duì)于蒙特卡洛采樣方法，當(dāng)采樣數(shù)達(dá)到2 000時(shí)才收斂，而3階Hermite多項(xiàng)式需要313個(gè)樣本就可以得到與蒙特卡洛采樣一致的結(jié)果?？紤]到蒙特卡洛采樣需要大量樣本，Liu等［13］利用代理模型技術(shù)，首先用100個(gè)訓(xùn)練樣本構(gòu)建Kriging代理模型，該模型可以代替耗時(shí)的CFD仿真預(yù)測(cè)RAE5243翼型上不同壓敏涂層厚度下的升力系數(shù)，誤差不超過0.6%，然后使用容易計(jì)算的代理模型進(jìn)行了107次蒙特卡洛采樣，獲得了幾何不確定性對(duì)翼型表面流動(dòng)的影響。相比而言，蒙特卡洛采樣直接模擬實(shí)際隨機(jī)過程，獲得的統(tǒng)計(jì)量精度高，可信度好，但樣本量需求十分龐大；混沌多項(xiàng)式展開將翼型性能參數(shù)的隨機(jī)響應(yīng)在混沌多項(xiàng)式基底上展開，利用多項(xiàng)式系數(shù)計(jì)算均值和標(biāo)準(zhǔn)差，樣本數(shù)需求有所降低，而精度與蒙特卡洛采樣相近。盡管如此，通過概率配點(diǎn)法求解多項(xiàng)式系數(shù)要求樣本數(shù)量不少于多項(xiàng)式基底數(shù)，限制了其在考慮高維不確定性因素的復(fù)雜系統(tǒng)中的實(shí)際應(yīng)用。

近年來(lái)，壓縮感知技術(shù)在信號(hào)處理領(lǐng)域取得迅速發(fā)展［14］，其利用原始信號(hào)的稀疏性質(zhì)，只需要很少的觀測(cè)樣本就可以用特定的重構(gòu)算法精確重構(gòu)原始信號(hào)［15］。針對(duì)不確定性量化中樣本數(shù)需求較多的問題，研究人員開始探索將壓縮感知技術(shù)應(yīng)用于不確定性量化中。Liang等［16］詳細(xì)探討了壓縮感知理論對(duì)研究對(duì)象的基本要求，包括系數(shù)的稀疏性和觀測(cè)矩陣的受限等距性質(zhì)（Restricted Isometry Property，RIP），指出混沌多項(xiàng)式展開在實(shí)際不確定性量化研究中滿足上述條件，并通過測(cè)試函數(shù)驗(yàn)證。Yang等［17］利用壓縮感知重構(gòu)混沌多項(xiàng)式系數(shù)，并在若干測(cè)試函數(shù)和偏微分方程組上測(cè)試，輸入變量維度涵蓋20～500范圍，所提出的方法即使在小樣本集上也能夠獲得良好的精度。陳江濤等［18］利用稀疏混沌多項(xiàng)式展開方法，在翼型繞流湍流模型稀疏不確定性和燒蝕材料物性參數(shù)不確定性2個(gè)算例上進(jìn)行研究，對(duì)于輸出量的統(tǒng)計(jì)特性稀疏混沌多項(xiàng)式展開能夠進(jìn)行準(zhǔn)確預(yù)測(cè)。Lu等［19］將壓縮感知中的正交匹配追蹤算法（Orthogonal Matching Pursuit，OMP）用于重構(gòu)混沌多項(xiàng)式系數(shù)，分別在考慮幾何不確定性的NLR7301翼型流動(dòng)和考慮壁面熱流不確定性的強(qiáng)迫對(duì)流換熱上對(duì)比驗(yàn)證，結(jié)果表明，OMP算法的優(yōu)勢(shì)在于達(dá)到相同的精度要求時(shí)收斂更快，觀測(cè)樣本量需求更少，適合復(fù)雜問題的不確定性分析。趙歡等［20］采用稀疏混沌多項(xiàng)式重構(gòu)算法，對(duì)高速自然層流翼型LRN1015開展考慮來(lái)流馬赫數(shù)和升力系數(shù)不確定性的穩(wěn)健性優(yōu)化，得益于自適應(yīng)前向-后向選擇的稀疏重構(gòu)算法，不確定性分析效率得到改進(jìn)，獲得了更好的穩(wěn)健性設(shè)計(jì)結(jié)果。盡管在壓縮感知技術(shù)與混沌多項(xiàng)式展開相結(jié)合方面國(guó)內(nèi)外學(xué)者已經(jīng)開展了初步探索，但在翼型幾何不確定性量化領(lǐng)域的研究仍相對(duì)較少。

本文考慮RAE2822翼型的幾何不確定性，將壓縮感知技術(shù)用于氣動(dòng)性能參數(shù)隨機(jī)響應(yīng)的不確定性量化，利用正交匹配追蹤算法重構(gòu)混沌多項(xiàng)式系數(shù)，并且將獲得的統(tǒng)計(jì)特性與傳統(tǒng)的蒙特卡洛采樣和滿秩概率配點(diǎn)法進(jìn)行對(duì)比分析，探討正交匹配追蹤重構(gòu)算法在收斂性能和樣本需求等方面的優(yōu)勢(shì)，為壓縮感知技術(shù)應(yīng)用于翼型不確定性量化領(lǐng)域的進(jìn)一步研究提供參考。

1 不確定性量化方法

1.1 混沌多項(xiàng)式展開

不確定性量化中的混沌多項(xiàng)式展開（Polynomial Chaos Expansion，PCE）方法將系統(tǒng)不確定性輸入和隨機(jī)響應(yīng)之間的映射關(guān)系在一系列正交多項(xiàng)式基底上展開。盡管現(xiàn)實(shí)中大多數(shù)信號(hào)都不具備稀疏性質(zhì)，但經(jīng)過混沌多項(xiàng)式展開后的系數(shù)將逐漸衰減，因此可以認(rèn)為具有一定的稀疏性。

混沌多項(xiàng)式展開可以表示為［21］

式中：y為系統(tǒng)隨機(jī)響應(yīng)輸出；x為系統(tǒng)不確定性輸入；Hi為與不確定性分布相對(duì)應(yīng)的混沌多項(xiàng)式基底。在本文研究范圍內(nèi)，假定幾何不確定性因素滿足正態(tài)分布，對(duì)應(yīng)的多項(xiàng)式基底為Hermite基底，其形式可以按照遞推公式［22］推導(dǎo)；ai為第i項(xiàng)基底對(duì)應(yīng)的多項(xiàng)式系數(shù)，需要通過采樣獲得若干組系統(tǒng)輸入-輸出響應(yīng)值后，利用最小二乘法求解：

其中：H可以寫為矩陣形式，行數(shù)為采樣樣本數(shù)N，列數(shù)（Npc+1）為總共的多項(xiàng)式項(xiàng)數(shù)：

Npc與不確定性輸入變量的維數(shù)m和多項(xiàng)式基底階數(shù)p有關(guān)：

從式（4）可見，隨著輸入維數(shù)的多項(xiàng)式基底階數(shù)的增加，基底總數(shù)呈階乘增長(zhǎng)，當(dāng)不確定性因素較多、系統(tǒng)響應(yīng)非線性較強(qiáng)時(shí)，多項(xiàng)式基底個(gè)數(shù)劇增，求解系數(shù)所需的樣本數(shù)也隨之增多，導(dǎo)致計(jì)算成本過高，即面臨“維數(shù)災(zāi)難”的困境。

1.2 壓縮感知理論

在求解多項(xiàng)式系數(shù)a時(shí)，研究人員希望需要的樣本量盡可能少，以提高不確定性量化效率。當(dāng)樣本數(shù)N小于多項(xiàng)式基底項(xiàng)數(shù)時(shí)，式（3）中的H矩陣行數(shù)小于列數(shù)，式（1）的求解轉(zhuǎn)變?yōu)榍范ǚ匠探M問題，壓縮感知理論正是關(guān)注此類問題的求解。從線性方程組的角度考慮，欠定方程組沒有唯一解，但壓縮感知理論證明，如果待求解向量滿足一定的條件，例如具有稀疏性，那么通過特定的重構(gòu)算法可以獲得唯一解。對(duì)于不確定性量化問題而言，如果系統(tǒng)響應(yīng)在混沌多項(xiàng)式基底上展開的系數(shù)足夠稀疏，那么系數(shù)向量可以通過如下l0范數(shù)最小化問題來(lái)求解：

式中：H矩陣起到觀測(cè)矩陣的作用；||a||0表示系數(shù)向量a中非零元素的個(gè)數(shù)，反映了系數(shù)向量的稀疏程度；ε為允許的擬合容差。

然而，直接求解上述l0范數(shù)最小化問題是一個(gè)NP-hard難題。為此，開發(fā)了一系列重構(gòu)算法，以可接受的計(jì)算成本獲得近似的最稀疏解，包括凸優(yōu)化算法、貪婪類算法等。其中，正交匹配追蹤（Orthogonal Matching Pursuit，OMP）是貪婪類重構(gòu)算法的代表，其主要優(yōu)勢(shì)是計(jì)算成本低、收斂速度快，適合交叉驗(yàn)證。

1.3 正交匹配追蹤算法

重構(gòu)算法是壓縮感知能否精確重構(gòu)原始信號(hào)的核心，根據(jù)應(yīng)用場(chǎng)景選擇合適的重構(gòu)算法至關(guān)重要。本文采用貪婪重構(gòu)算法中的正交匹配追蹤重構(gòu)混沌多項(xiàng)式系數(shù)。

由于遍歷H矩陣中所有列向量基底組合過于復(fù)雜，正交匹配追蹤算法采用貪婪思想，每次從列向量基底中選擇與殘差最相關(guān)的基底，添加到已有的支持集中，在選出的基底上嘗試重構(gòu)多項(xiàng)式系數(shù)，系數(shù)向量中只有選出的基底對(duì)應(yīng)的系數(shù)不為零，而其余系數(shù)均為零，以此實(shí)現(xiàn)重構(gòu)系數(shù)的稀疏性。上述步驟迭代若干次，直至殘差降至指定的容差范圍內(nèi)，算法流程如圖1所示。

圖1 正交匹配追蹤算法流程Fig.1 Flowchart of OMP

由于在每一個(gè)迭代步，OMP算法都是尋找能使殘差下降最快的基向量，因此能夠以較快的速度收斂到唯一的稀疏解。

通過重構(gòu)算法可以獲得系數(shù)向量a，進(jìn)而獲得混沌多項(xiàng)式展開完整的表達(dá)式，可以根據(jù)式（1）計(jì)算系統(tǒng)輸出響應(yīng)。對(duì)于不確定性量化問題，研究人員更加關(guān)注隨機(jī)響應(yīng)的統(tǒng)計(jì)特性，包括均值μ和標(biāo)準(zhǔn)差σ，兩者可以直接通過混沌多項(xiàng)式系數(shù)分量來(lái)估計(jì)：

1.4 標(biāo)準(zhǔn)測(cè)試算例驗(yàn)證

綜合上述方法，建立了一個(gè)可以減少樣本數(shù)需求的不確定性量化框架，如圖2所示。首先確定不確定性的輸入維度和概率分布，然后根據(jù)概率分布和多項(xiàng)式階數(shù)生成混沌多項(xiàng)式基底，所需的樣本從高斯積分點(diǎn)中隨機(jī)選取，確定觀測(cè)矩陣H。對(duì)于翼型流動(dòng)問題，通過計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）仿真獲得各個(gè)樣本下的氣動(dòng)力系數(shù)和其他關(guān)注的流場(chǎng)量，最后在建立的樣本集上使用正交匹配追蹤重構(gòu)算法，求解混沌多項(xiàng)式系數(shù)，并獲得所關(guān)注流動(dòng)參數(shù)的統(tǒng)計(jì)特性。

圖2 不確定性量化流程Fig.2 Uncertainty quantification framework

為了驗(yàn)證上述不確定性量化框架的精度，本文首先在Ishigami測(cè)試函數(shù)上進(jìn)行算例驗(yàn)證。該函數(shù)具有很強(qiáng)的非線性，在不確定性量化研究中廣泛用于精度校驗(yàn)［23］。

該函數(shù)的定義為

其均值μ和標(biāo)準(zhǔn)差σ可以通過理論分析精確計(jì)算，作為評(píng)估不確定度量化誤差的基值。將上述不確定性量化框架用于該測(cè)試算例，并與理論值和傳統(tǒng)的不確定性量化方法，包括蒙特卡洛采樣和滿秩概率配點(diǎn)法（Full Rank Probabilistic Collocation，F(xiàn)RPC）進(jìn)行對(duì)比。

圖3給出了均值和標(biāo)準(zhǔn)差相對(duì)誤差隨樣本數(shù)增加的變化過程。相比之下，3種方法中蒙特卡洛采樣收斂最慢，重構(gòu)誤差隨著樣本數(shù)的增加緩慢降低?；煦缍囗?xiàng)式展開基底輸入維數(shù)為3，多項(xiàng)式階數(shù)取13，根據(jù)式（4）可知共有560個(gè)基底，當(dāng)變量數(shù)小于基底數(shù)時(shí)，方程組欠定，難以獲得精確解，故均值和標(biāo)準(zhǔn)差的重構(gòu)誤差都較大；而當(dāng)樣本數(shù)等于基底數(shù)時(shí)，方程組滿秩，可以用最小二乘法求得精確解，因而重構(gòu)誤差曲線陡降。而采用同樣混沌多項(xiàng)式基底的正交匹配追蹤算法，僅需140個(gè)樣本即可收斂到與滿秩概率配點(diǎn)法相同的精度水平，這是因?yàn)?60個(gè)混沌多項(xiàng)式基底對(duì)應(yīng)的系數(shù)向量具有稀疏性，即大部分系數(shù)幅值接近于0，表明這些系數(shù)對(duì)應(yīng)的多項(xiàng)式基底影響作用較小，而正交匹配追蹤算法可以優(yōu)先篩選占主要影響的基底并求解對(duì)應(yīng)的系數(shù)，從而實(shí)現(xiàn)用較少的樣本獲得精確的不確定性量化結(jié)果。

圖3 蒙特卡洛采樣、滿秩概率配點(diǎn)和正交匹配追蹤算法對(duì)Ishigami測(cè)試算例的預(yù)測(cè)誤差Fig.3 Prediction errors of Ishigami function by MCS，F(xiàn)RPC and OMP

2 RAE2822不確定性量化分析

在通過測(cè)試算例驗(yàn)證所用方法的精度后，將上述框架用于RAE2822翼型的幾何不確定性量化分析。首先對(duì)所采用的數(shù)值模型進(jìn)行驗(yàn)證，將仿真結(jié)果與試驗(yàn)結(jié)果對(duì)比。為了描述翼型幾何不確定性變化，采用類函數(shù)/形函數(shù)變換（Class-Shape Transformation，CST）將 原 始 翼 型 參 數(shù)化，將翼型幾何的不確定性變化轉(zhuǎn)變?yōu)镃ST參數(shù)的隨機(jī)波動(dòng)，并確定各參數(shù)概率分布規(guī)律和變化范圍。最后，在本文所研究的翼型幾何不確定性因素下，用混沌多項(xiàng)式展開和正交匹配追蹤算法開展不確定性量化，獲得翼型升力系數(shù)、阻力系數(shù)和其他主要流場(chǎng)量的統(tǒng)計(jì)特性，并與蒙特卡洛采樣和概率配點(diǎn)法在收斂性、量化精度和樣本數(shù)需求等方面進(jìn)行對(duì)比，以驗(yàn)證所用方法的有效性。

2.1 數(shù)值模型驗(yàn)證

RAE2822翼型繞流是典型的二維跨聲速可壓縮流動(dòng)，被歐洲EUROVAL項(xiàng)目和AGARD選為經(jīng)典驗(yàn)證算例。Cook等［24］在RAE 8 ft×6 ft（1 ft=0.304 8 m）跨聲速風(fēng)洞中進(jìn)行流場(chǎng)實(shí)驗(yàn)測(cè)量，共選擇10個(gè)工況，通過壓力探針測(cè)量壓力分布，而升力、阻力和力矩系數(shù)通過對(duì)翼型表面壓力進(jìn)行積分獲得。此外，還用邊界層探針測(cè)量了邊界層速度剖面。

RAE2822翼型的幾何形狀以上下表面離散點(diǎn)的形式表示，為了描述翼型幾何形狀的不確定性變化，需要用參數(shù)化方法將離散的翼型幾何點(diǎn)轉(zhuǎn)化為控制參數(shù)。本文采用Kulfan［25］提出的類形狀/形函數(shù)變換方法來(lái)擬合原始翼型。CST方法使用類函數(shù)C（φ）和形函數(shù)S（φ）來(lái)定義翼型的幾何形狀：

式中：φ和ζ是基于弦長(zhǎng)歸一化的翼型坐標(biāo)；ζT是尾緣厚度，對(duì)于RAE2822類型的圓頭尖尾翼型，尾緣厚度為0。

C（φ）是用于定義翼型整體類型的類函數(shù)，其表達(dá)式為

通過調(diào)整N1和N2，可以生成不同類型的翼型，如圖4所示。

圖4 類函數(shù)C（φ）決定的不同N1、N2取值對(duì)應(yīng)的基本翼型Fig.4 Basic airfoil determined by class function C（φ）with different N1， N2

對(duì)于前緣圓形、尾緣尖形的RAE2822翼型，N1取0.5，N2取1。形函數(shù)S（φ）通過多階Bernstein多項(xiàng)式的疊加來(lái)描述翼型的具體形狀：

式中：Ai是第i項(xiàng)的權(quán)重系數(shù)；n是Bernstein多項(xiàng)式 階 數(shù)（Bernstein Polynomial Order，BPO）。CST參數(shù)化方法即是通過控制形函數(shù)中的Ai來(lái)擬合成不同形狀的翼型，擬合誤差隨著多項(xiàng)式階數(shù)的增加而減小，但與此同時(shí)變量Ai的個(gè)數(shù)也隨之增加。

圖5給出了不同階數(shù)下的對(duì)原始翼型的擬合誤差，可見4階CST已經(jīng)能夠在風(fēng)洞模型公差范圍內(nèi)擬合RAE2822翼型，可以滿足實(shí)際應(yīng)用需要。

圖5 不同Bernstein基函數(shù)階數(shù)的CST方法對(duì)RAE2822翼型的擬合誤差Fig.5 Fitting error of RAE2822 airfoil using CST with different BPO

采用4階CST參數(shù)化，上下表面各需要5個(gè)參數(shù)（A0～A9）來(lái)描述無(wú)量綱翼型形狀；此外，為了保持與實(shí)驗(yàn)?zāi)Ｐ鸵恢乱则?yàn)證數(shù)值模型精度，同時(shí)考慮弦長(zhǎng)不確定性的變化，翼型的弦長(zhǎng)c也作為幾何參數(shù)。因此，總共確定了11個(gè)幾何參數(shù)。

為了驗(yàn)證數(shù)值模型，首先進(jìn)行確定性仿真。在不確定性量化研究中，通常采用來(lái)流馬赫數(shù)Ma為0.734，攻角為2.79°，對(duì)應(yīng)雷諾數(shù)Re為6.5×106的工況［26］。在數(shù)值仿真驗(yàn)證時(shí)，選取與其最為接近的實(shí)驗(yàn)測(cè)量工況進(jìn)行仿真，對(duì)應(yīng)的Ma為0.730，攻角為3.19°。為了同時(shí)驗(yàn)證經(jīng)過CST擬合的翼型氣動(dòng)性能，對(duì)原始翼型和擬合翼型均進(jìn)行雷諾平均納維-斯托克斯方程（Reynold-Averaged Navier-Stokes equations，RANS）模擬，計(jì)算域邊界距離翼型表面25倍弦長(zhǎng)，以避免計(jì)算域邊界對(duì)翼型附近流動(dòng)的干擾。對(duì)計(jì)算域劃分C型拓?fù)渚W(wǎng)格，邊界層第1層網(wǎng)格高度由y+值控制，本研究采用SSTk-ω湍流模型，要求壁面上y+＜1。網(wǎng)格節(jié)點(diǎn)數(shù)通過網(wǎng)格無(wú)關(guān)性驗(yàn)證確定，進(jìn)行了一系列網(wǎng)格單元總數(shù)逐漸遞增的仿真，表1列出了升力系數(shù)（Lift coefficient，CL）和阻力系數(shù)（Drag coefficient，CD）與實(shí)驗(yàn)結(jié)果的對(duì)比。表中結(jié)果表明，當(dāng)網(wǎng)格單元數(shù)達(dá)到120 000后，計(jì)算結(jié)果隨網(wǎng)格單元的進(jìn)一步增加變化不大。因此，在后續(xù)仿真將基于129 258單元數(shù)的網(wǎng)格開展，以平衡仿真精度和時(shí)間成本。該套網(wǎng)格在翼型附近區(qū)域的分布情況如圖6所示，在翼型上表面中間弦長(zhǎng)后的激波存在區(qū)域進(jìn)行網(wǎng)格節(jié)點(diǎn)加密，以提高對(duì)激波的辨識(shí)度。

表1 網(wǎng)格無(wú)關(guān)性驗(yàn)證Table 1 Grid independence test

圖6 RAE2822翼型計(jì)算網(wǎng)格Fig.6 Computation grid around RAE2822 airfoil

在確定了網(wǎng)格單元數(shù)后，用相同的網(wǎng)格拓?fù)鋵?duì)擬合翼型也劃分了計(jì)算域網(wǎng)格，圖7中給出了原始翼型和擬合翼型表面壓力系數(shù)（Pressure coefficient，Cp）分布與實(shí)驗(yàn)結(jié)果的對(duì)比，結(jié)果表明，經(jīng)過CST擬合的翼型仿真結(jié)果與原始翼型取得一致，從仿真結(jié)果層面驗(yàn)證了擬合精度。

圖7 原始翼型和擬合翼型壓力系數(shù)的對(duì)比Fig.7 Comparison of original and fitting airfoils on pressure coefficient

2.2 幾何不確定性范圍和概率分布

在驗(yàn)證數(shù)值模型的基礎(chǔ)上，可以開展考慮幾何不確定性的RAE2822翼型氣動(dòng)性能不確定性量化。由于翼型幾何已經(jīng)通過CST方法進(jìn)行參數(shù)化，因此翼型幾何不確定性的變化可以轉(zhuǎn)化為CST參數(shù)的隨機(jī)波動(dòng)。參考其他有關(guān)翼型幾何不確定性的研究［27］，針對(duì)加工、裝配和實(shí)際飛行工況等因素導(dǎo)致的實(shí)際翼型與原始翼型差異，將翼型弦長(zhǎng)的不確定性變化也考慮在內(nèi)，確定了本研究中RAE2822的幾何厚度變化范圍為±1 mm。圖8中給出翼型前緣和尾緣的不確定性變化邊界。

圖8 RAE2822在前緣和尾緣的幾何不確定性變化范圍Fig.8 Geometrical uncertainty range of RAE2822 airfoil at leading edge and trailing edge

11個(gè)CST幾何參數(shù)的范圍由圖中最薄和最厚的翼型確定，這些參數(shù)假設(shè)均服從高斯分布N（μ，σ）。為了實(shí)現(xiàn)CST參數(shù)隨機(jī)變化時(shí)翼型幾何能夠保持在最薄和最厚翼型之間，采用最小二乘法分別反求加厚翼型和減薄翼型對(duì)應(yīng)的CST參數(shù)值，由此確定了11個(gè)CST參數(shù)能夠變化的上下限范圍。根據(jù)3σ原則，假定此范圍為±3σ的邊界，由此確定了各變量標(biāo)準(zhǔn)差σ的取值，均值μ則由原始翼型反求CST參數(shù)獲得。各CST參數(shù)的均值μ和標(biāo)準(zhǔn)差σ在表2中列出。

表2 11個(gè)CST參數(shù)的均值μ和標(biāo)準(zhǔn)差σTable 2 Mean μ and standard deviation σ of 11 CST parameters

圖9給出了按照上述CST參數(shù)化方法生成的一系列隨機(jī)翼型，這些翼型分布在加厚翼型和減薄翼型之間，其分布規(guī)律表明所采用的方法生成的隨機(jī)翼型可以覆蓋±1 mm的幾何不確定性變化區(qū)間，可以用于描述翼型幾何不確定性變化。

圖9 隨機(jī)翼型幾何分布圖Fig.9 Distribution of random airfoils geometry

幾何不確定性變化范圍確定后，采用圖2所示的不確定性量化框架，用正交匹配追蹤算法求解混沌多項(xiàng)式系數(shù)并獲得RAE2822翼型氣動(dòng)性能的統(tǒng)計(jì)特性，與蒙特卡洛采樣（MCS）和滿秩概率配點(diǎn)法（FRPC）對(duì)比。

2.3 氣動(dòng)力系數(shù)對(duì)比

表3中列出了3種方法求解的升力系數(shù)和阻力系數(shù)的統(tǒng)計(jì)特性。由于對(duì)于實(shí)際翼型的不確定性量化問題尚無(wú)理論解，考慮到蒙特卡洛采樣法模擬真實(shí)概率分布采樣，樣本信息豐富，因此以蒙特卡洛采樣獲得的結(jié)果為參考計(jì)算相對(duì)誤差。

表3 蒙特卡洛采樣、滿秩概率配點(diǎn)法和正交匹配追蹤預(yù)測(cè)的氣動(dòng)力系數(shù)均值μ和標(biāo)準(zhǔn)差σTable 3 Mean μ and standard deviation σ of aerodynamic coefficients predicted by MCS， FRPC and OMP

由表3可知，在幾何不確定性作用下，翼型升力系數(shù)均值相比于原型下降了約1%，表明幾何不確定性的存在導(dǎo)致翼型升力降低；另一方面，幾何不確定性導(dǎo)致升力系數(shù)和阻力系數(shù)出現(xiàn)一定程度的波動(dòng)，其中阻力系數(shù)變化的標(biāo)準(zhǔn)差可以達(dá)到均值的4%。比較而言，滿秩概率配點(diǎn)法所獲得的升力系數(shù)和阻力系數(shù)與蒙特卡洛采樣較為接近，但需要364個(gè)樣本才能使H矩陣滿秩；而正交匹配追蹤算法只需要120個(gè)樣本，即可取得與滿秩概率配點(diǎn)法基本一致的結(jié)果，表明壓縮感知方法可以用相對(duì)較少的樣本獲得較高的精度。

上述結(jié)果體現(xiàn)了樣本量充足時(shí)各方法的預(yù)測(cè)精度，而圖10中則給出不同樣本數(shù)下的相對(duì)誤差來(lái)比較各方法的收斂性，其中相對(duì)誤差是以最多樣本數(shù)下蒙特卡洛采樣的結(jié)果為基準(zhǔn)計(jì)算的。3種方法中，蒙特卡洛采樣收斂速度最慢，幾乎在所有樣本數(shù)水平下相對(duì)誤差都是最高的?；诟怕逝潼c(diǎn)法的PCE方法直到樣本數(shù)達(dá)到364，即滿足滿秩條件時(shí)才快速收斂，364個(gè)樣本的結(jié)果與1 200個(gè)樣本量的蒙特卡洛采樣方法相比，相對(duì)誤差減少了一個(gè)數(shù)量級(jí)；而在樣本數(shù)小于多項(xiàng)式系數(shù)個(gè)數(shù)時(shí)，由于方程組欠定，求解誤差較大。

圖10 蒙特卡洛采樣、滿秩概率配點(diǎn)和正交匹配追蹤算法收斂情況對(duì)比Fig.10 Convergence comparison of MCS， FRPC and OMP

比較而言，正交匹配追蹤算法可以用更少的樣本獲得與概率配點(diǎn)法相接近的精度，特別是當(dāng)樣本數(shù)低于滿秩要求，相比于概率配點(diǎn)法用最小二乘擬合多項(xiàng)式系數(shù)，正交匹配追蹤的誤差更低或持平。不過，在樣本數(shù)增加到一定規(guī)模之前，正交匹配追蹤算法和概率配點(diǎn)法之間的差距并不明顯。具體對(duì)于本算例而言，均值μ的預(yù)測(cè)出現(xiàn)顯著差異是在樣本數(shù)達(dá)到30以上，而對(duì)于標(biāo)準(zhǔn)差σ則是60個(gè)樣本，表明標(biāo)準(zhǔn)差的準(zhǔn)確預(yù)測(cè)比均值難度更大。

正交匹配追蹤算法能夠表現(xiàn)出上述優(yōu)勢(shì)得益于多項(xiàng)式系數(shù)的稀疏性。盡管在實(shí)際應(yīng)用中假設(shè)了稀疏性的前提，但重構(gòu)獲得的多項(xiàng)式系數(shù)是否足夠稀疏仍然需要校驗(yàn)。圖11對(duì)比了2種PCE方法中多項(xiàng)式基底對(duì)應(yīng)系數(shù)的幅值，按基底序號(hào)排序。圖中的2組系數(shù)分別是364個(gè)樣本的滿秩概率配點(diǎn)法PCE和120個(gè)樣本的正交匹配追蹤獲得的，為了直觀顯示多項(xiàng)式系數(shù)的相對(duì)大小，橫縱坐標(biāo)軸都采用對(duì)數(shù)刻度。圖中系數(shù)幅值分布表明，大部分PCE系數(shù)都接近于0，如果假定小于10-6的系數(shù)可被視為0，那么系數(shù)向量確實(shí)滿足稀疏性條件。從另一角度分析，所有通過正交匹配追蹤求解出的系數(shù)幅值都大于10-6，這表明該算法能夠從所有的多項(xiàng)式基底中篩選出起主導(dǎo)作用的子集，從而利用系數(shù)向量的稀疏性精確重構(gòu)對(duì)應(yīng)的非零系數(shù)。

圖11 滿秩概率配點(diǎn)和正交匹配追蹤重構(gòu)的混沌多項(xiàng)式系數(shù)對(duì)比Fig.11 Comparison of polynomial chaos coefficients reconstructed by FRPC and OMP

除了對(duì)比均值和標(biāo)準(zhǔn)差等統(tǒng)計(jì)特性之外，還需要進(jìn)一步驗(yàn)證概率密度分布。圖12中給出了采用蒙特卡洛采樣、滿秩法PCE和正交匹配追蹤預(yù)測(cè)的升力系數(shù)和阻力系數(shù)的概率密度分布函數(shù)（Probability Density Function，PDF）。圖中結(jié)果表明，滿秩概率配點(diǎn)法PCE和正交匹配追蹤能夠精確重構(gòu)出與蒙特卡洛采樣結(jié)果一致的氣動(dòng)力系數(shù)概率分布，且正交匹配追蹤只需要120個(gè)觀測(cè)樣本，遠(yuǎn)少于需要364個(gè)樣本才能滿秩的PCE和3 000個(gè)樣本才能收斂的蒙特卡洛采樣。

圖12 蒙特卡洛采樣、滿秩概率配點(diǎn)和正交匹配追蹤算法預(yù)測(cè)的升力系數(shù)、阻力系數(shù)概率密度函數(shù)對(duì)比Fig.12 Comparison of probability density function of CL and CD predicted by MCS， FRPC and OMP

綜上，在考慮幾何不確定性的翼型氣動(dòng)力系數(shù)的預(yù)測(cè)中，壓縮感知技術(shù)在樣本需求量少方面展現(xiàn)出顯著優(yōu)勢(shì)，能夠高效靈活地進(jìn)行不確定性量化分析。

2.4 流場(chǎng)參數(shù)對(duì)比

氣動(dòng)力系數(shù)是將翼型表面所受壓力積分獲得的，為了進(jìn)一步對(duì)比流場(chǎng)重構(gòu)結(jié)果，提取翼型上下表面壓力系數(shù)沿弦長(zhǎng)的分布，如圖13所示。

圖13 翼型表面壓力系數(shù)沿弦長(zhǎng)的均值和標(biāo)準(zhǔn)差分布Fig.13 Mean and standard deviation of Cp on airfoil surfaces along chord

RAE2822作為典型的超臨界翼型，氣流在翼型上表面加速降壓，從前緣起有一個(gè)較為穩(wěn)定的低壓區(qū)，直到出現(xiàn)激波后壓力驟升，而在下表面上壓力系數(shù)沒有突變，且在翼型尾緣附近壓力有所升高，以提供更大升力。在均值方面，滿秩概率配點(diǎn)法和正交匹配追蹤都可以獲得與蒙特卡洛采樣相同的壓力系數(shù)分布規(guī)律，結(jié)果一致性較好；而在標(biāo)準(zhǔn)差方面，對(duì)于壓力波動(dòng)不明顯的區(qū)域，包括翼型上表面沒有激波的區(qū)域和整個(gè)翼型下表面，標(biāo)準(zhǔn)差σ的預(yù)測(cè)結(jié)果也很接近，3種方法的主要差異體現(xiàn)在60%弦長(zhǎng)附近的激波所在區(qū)域。幾何不確定性的存在使得翼型偏離設(shè)計(jì)形狀，不但會(huì)使激波的位置發(fā)生改變，激波強(qiáng)度也會(huì)因波前氣流加速程度不同而變化。由于氣流經(jīng)過激波后壓力驟升，激波強(qiáng)度和位置的變化將導(dǎo)致激波所在區(qū)域壓力系數(shù)變化劇烈，即標(biāo)準(zhǔn)差增大，同時(shí)不同不確定性量化方法之間標(biāo)準(zhǔn)差預(yù)測(cè)的差異也被放大。相比而言，在標(biāo)準(zhǔn)差最大的位置，正交匹配追蹤的結(jié)果更加接近蒙特卡洛采樣的參考值，而滿秩法PCE的結(jié)果則高估了壓力系數(shù)的波動(dòng)程度。除了激波所在區(qū)域之外，翼型前緣壓力波動(dòng)也較為顯著，這是因?yàn)橐硇拖卤砻骜v點(diǎn)位置因幾何形狀發(fā)生微小改變而移動(dòng)。

幾何不確定性的存在不僅影響著翼型表面壓力，同時(shí)也改變了翼型附近的其他流場(chǎng)參數(shù)。在圖14中，給出了蒙特卡洛采樣、滿秩法和正交匹配追蹤算法重構(gòu)的馬赫數(shù)均值和標(biāo)準(zhǔn)差分布，與前述的氣動(dòng)力系數(shù)不同，流場(chǎng)馬赫數(shù)的重構(gòu)是在每個(gè)網(wǎng)格節(jié)點(diǎn)上進(jìn)行的。正交匹配追蹤算法對(duì)馬赫數(shù)均值的預(yù)測(cè)在大多數(shù)位置上與蒙特卡洛采樣非常接近，上表面空氣從前緣開始加速至超聲速，經(jīng)過激波后變?yōu)閬喡曀伲R赫數(shù)驟變。標(biāo)準(zhǔn)差的差異主要體現(xiàn)在60%弦長(zhǎng)附近的激波區(qū)域，這也與壓力系數(shù)分布呈現(xiàn)的規(guī)律一致，幾何不確定性影響了激波出現(xiàn)的位置以及激波前后馬赫數(shù)差值，重構(gòu)激波區(qū)域內(nèi)馬赫數(shù)的概率特性也更具挑戰(zhàn)。正交匹配追蹤的結(jié)果表明，翼型附近區(qū)域的流場(chǎng)參數(shù)變化同樣可以通過稀疏重構(gòu)的方式以一定的精度準(zhǔn)確預(yù)測(cè)，且所需樣本量更少，從而在不確定性量化中展現(xiàn)出顯著的優(yōu)勢(shì)。

3 結(jié) 論

本文將壓縮感知技術(shù)與不確定性量化相結(jié)合，對(duì)正交匹配追蹤算法用于RAE2822翼型幾何不確定性量化進(jìn)行了研究，可以得到以下結(jié)論：

1）飛行器機(jī)翼在加工、裝配和實(shí)際飛行過程中，會(huì)因?yàn)楦鞣N因素存在幾何外形的不確定性，將RAE2822翼型經(jīng)過CST方法參數(shù)化后，翼型幾何的不確定性可以轉(zhuǎn)換為CST參數(shù)的隨機(jī)波動(dòng)。

2）經(jīng)過不確定性量化研究，幾何外形不確定性的存在會(huì)導(dǎo)致翼型升力系數(shù)降低，且翼型上表面激波的強(qiáng)度和位置發(fā)生變化，激波所在區(qū)域的流場(chǎng)參數(shù)出現(xiàn)波動(dòng)，影響飛行器的正常飛行。

3）通過測(cè)試函數(shù)和RAE2822翼型驗(yàn)證，與傳統(tǒng)的蒙特卡洛采樣、滿秩概率配點(diǎn)法等不確定性量化方法相比，正交匹配追蹤可以用相對(duì)較少的樣本獲得與相同精度水平的量化結(jié)果。

4）正交匹配追蹤算法用稀疏重構(gòu)的方式求解混沌多項(xiàng)式系數(shù)，不僅能夠預(yù)測(cè)氣動(dòng)力系數(shù)的統(tǒng)計(jì)特性和概率分布，還可以重構(gòu)翼型附近的流場(chǎng)參數(shù)變化規(guī)律，有助于進(jìn)一步研究不確定性對(duì)翼型流動(dòng)的影響。

通過測(cè)試函數(shù)和RAE2822翼型幾何不確定性的研究，驗(yàn)證了壓縮感知技術(shù)用于不確定性量化在效率、收斂性和精度等方面的優(yōu)勢(shì)。后續(xù)的工作將在此基礎(chǔ)上開展參數(shù)敏感度分析，進(jìn)一步研究不同幾何不確定性參數(shù)對(duì)翼型氣動(dòng)性能的影響程度，為開展翼型穩(wěn)健性設(shè)計(jì)提供參考。