李春娜，賈續(xù)毅，龔春林

（西北工業(yè)大學(xué)空天飛行技術(shù)研究所，西安710072）

0 引言

在氣動(dòng)優(yōu)化設(shè)計(jì)中，CFD方法由于具有高可信度而逐漸得到廣泛應(yīng)用［1-2］，但是大量的CFD分析所消耗的時(shí)間和計(jì)算資源也是十分巨大的。為了有效地降低優(yōu)化成本，研究人員提出了一系列更高效的優(yōu)化設(shè)計(jì)方法，例如，使用伴隨方法可以有效地提高多設(shè)計(jì)變量梯度優(yōu)化過程的效率［3-4］。但是從優(yōu)化算法本身進(jìn)行改進(jìn)仍不能有效地解決調(diào)用大量CFD的問題。

代理優(yōu)化算法［5-6］的應(yīng)用使CFD的調(diào)用次數(shù)呈量級(jí)式下降，大幅提高了優(yōu)化效率。例如，韓少?gòu)?qiáng)等［7］使用梯度增強(qiáng)型Kriging模型實(shí)現(xiàn)了高維設(shè)計(jì)變量的氣動(dòng)反設(shè)計(jì)；邱亞松等［8］使用本征正交分解（Proper Orthogonal Decomposition，簡(jiǎn)稱POD）和Kriging/RBF代理模型實(shí)現(xiàn)了對(duì)變外形的翼型定常流場(chǎng)的預(yù)測(cè)。但是，Kriging等代理模型一般只適用于單輸出，對(duì)于流場(chǎng)參數(shù)預(yù)測(cè)等多輸出問題建模效率低下。

近年來，研究人員將機(jī)器學(xué)習(xí)方法應(yīng)用于流場(chǎng)預(yù)測(cè)、氣動(dòng)優(yōu)化等方面，并取得顯著成效［9-11］。Zhu Linyang等［12］使用RBF神經(jīng)網(wǎng)絡(luò)建立了數(shù)據(jù)驅(qū)動(dòng)的湍流模型并實(shí)現(xiàn)了渦黏的預(yù)測(cè)；劉凌君等［13］使用反向傳播神經(jīng)網(wǎng)絡(luò)（Back Propagation based Neural Network，簡(jiǎn)稱BPNN）構(gòu)建了翼型參數(shù)化系數(shù)與翼面壓力系數(shù)的神經(jīng)網(wǎng)絡(luò)，并實(shí)現(xiàn)了翼型反設(shè)計(jì)，但存在模型訓(xùn)練樣本量大（4 000個(gè)訓(xùn)練樣本），需要將翼型參數(shù)化系數(shù)分組建立多個(gè)神經(jīng)網(wǎng)絡(luò)等復(fù)雜問題。

直接使用神經(jīng)網(wǎng)絡(luò)（Neural Network，簡(jiǎn)稱NN）構(gòu)建流場(chǎng)預(yù)測(cè)模型，由于流場(chǎng)數(shù)據(jù)維度高，模型輸入輸出維度差異過大，導(dǎo)致模型學(xué)習(xí)緩慢、建模困難。而以POD為代表的降階模型可以有效地解決高維數(shù)據(jù)的維度災(zāi)難問題，并已成功應(yīng)用于流場(chǎng)近似求解、翼型反設(shè)計(jì)等方面［14-15］。

本文以翼型為研究對(duì)象，首先，結(jié)合POD和NN的特點(diǎn)，發(fā)展一種基于兩層POD和BPNN的翼型反設(shè)計(jì)方法；然后，將該方法應(yīng)用于亞/跨聲速下的翼型反設(shè)計(jì)，包括翼型庫(kù)建立、聚類取樣、基模態(tài)個(gè)數(shù)選取、模型訓(xùn)練等方面；最后，通過算例對(duì)預(yù)測(cè)誤差、聚類取樣效果和超參數(shù)等方面進(jìn)行分析。

1 建模方法

1.1 樣本庫(kù)建立

常見的翼型參數(shù)化方法包括Hicks-Henne、CST、PARSEC、B樣條等［16-17］，其中Hicks-Henne參數(shù)化通過在基準(zhǔn)翼型的上下表面疊加數(shù)個(gè)Hicks-Henne形函數(shù)，并調(diào)節(jié)形函數(shù)的個(gè)數(shù)、作用位置、高度系數(shù)來實(shí)現(xiàn)翼型外形的改變，在翼型反設(shè)計(jì)、氣動(dòng)優(yōu)化設(shè)計(jì)中得到廣泛應(yīng)用。

本文以典型的亞聲速翼型NACA0012和跨聲速翼型RAE2822為基準(zhǔn)翼型，在其上下翼面各施加6個(gè)Hicks-Henne形函數(shù)，其高度系數(shù)范圍在亞聲 速 中 為［±0.005，±0.010，±0.010，±0.010，±0.010，±0.005］，在 跨 聲 速 中 為［±0.005，±0.005，±0.005，±0.005，±0.005，±0.005］，對(duì)應(yīng)的作用位置均為0.10c、0.25c、0.45c、0.65c、0.80c、0.90c（c為翼型弦長(zhǎng)）。

1.2 POD降階

POD是將具有高維特征的數(shù)據(jù)轉(zhuǎn)化為正交的低維特征數(shù)據(jù)。為了最大程度地保留原始數(shù)據(jù)的信息量，應(yīng)使降階后的數(shù)據(jù)方差最大化［18］。

設(shè)r個(gè)樣本構(gòu)成原始樣本數(shù)據(jù)矩陣X=每 個(gè) 樣 本 為n維 向 量，即對(duì)X做中心化處理，得到新樣本數(shù)據(jù)D=其中

使用線性映射矩陣U（m×n）將D降 階 至m維 空間，降階所得數(shù)據(jù)為Y（m×r），即

要使Y方差最大化，只需計(jì)算D的協(xié)方差矩陣var(D)的特征值和特征向量，并將前m個(gè)最大的特征值所對(duì)應(yīng)的特征向量q1，q2，…，qm構(gòu)成U即可。var(D)的計(jì)算公式為

式中：DT為D的轉(zhuǎn)置矩陣。

POD降價(jià)流程如圖1所示。

圖1 POD降階流程Fig.1 Process of building POD reduced-order model

1.3 BPNN模型

BPNN是一種通過前向傳播計(jì)算輸出，反向傳播計(jì)算誤差的全連接神經(jīng)網(wǎng)絡(luò)。其特點(diǎn)是多輸入多輸出建模，具備較強(qiáng)的非線性擬合能力，普遍應(yīng)用于分類或回歸問題［19］。典型的BPNN具有三層：輸入層、隱藏層和輸出層，如圖2所示。

圖2 三層BPNNFig.2 A BPNN with three-layer

對(duì)于如圖2所示的BPNN模型，輸入層為m個(gè)神經(jīng)元，輸出層為n個(gè)神經(jīng)元，隱藏層神經(jīng)元個(gè)數(shù)為q，其大小由m、n值確定，這里q=max{m，n}。每層神經(jīng)元通過權(quán)值ω和激活函數(shù)連接［19］，本文算例中隱藏層和輸出層的激活函數(shù)類型分別為tansig和purelin函 數(shù)，BPNN的 訓(xùn) 練 流 程 如 表1所示。

表1 BPNN的訓(xùn)練流程Table 1 Training process of BPNN

1.4 翼型反設(shè)計(jì)方法

結(jié)合POD和BPNN模型，發(fā)展一種快速的翼型反設(shè)計(jì)方法，其流程圖如圖3所示。具體流程分為四步。

圖3 基于兩層POD和BPNN的翼型反設(shè)計(jì)方法Fig.3 Airfoil reverse design method based on two-layer POD and BPNN

第一步，翼型庫(kù)建立。

使用參數(shù)化方法得到r個(gè)不同幾何外形的翼型（無量綱化），每個(gè)翼型可以用坐標(biāo)集表示，其中x0為構(gòu)成翼型的所有坐標(biāo)點(diǎn)的橫坐標(biāo)，每個(gè)翼型均相同；yi為構(gòu)成第i個(gè)翼型的所有坐標(biāo)點(diǎn)的縱坐標(biāo)。因此，每個(gè)翼型可以由唯一的yi確定。使用流場(chǎng)求解器，可以計(jì)算給定工況下每個(gè)翼型所對(duì)應(yīng)的翼面壓力系數(shù)。第i個(gè)翼型對(duì)應(yīng)的翼面壓力系數(shù)用{xC，Ci}表示，其中xC為描述上下翼面所有點(diǎn)的橫坐標(biāo)值，且每個(gè)翼型均一致；Ci為第i個(gè)翼型的翼面壓力系數(shù)值。因此，每個(gè)翼型的翼面壓力系數(shù)可以由唯一的Ci確定。

第二步，確定樣本集。

在建立翼面壓力系數(shù)到幾何外形的映射之前，為了確定最佳訓(xùn)練樣本量，提高訓(xùn)練效率，本文提出一種聚類取樣策略。首先將所有樣本的翼面壓力系數(shù)C進(jìn)行POD降階，得到壓力系數(shù)的基模態(tài)系數(shù)矩陣將基模態(tài)系數(shù)一致的樣 本 歸為一類Si(i=1，2，…，k，且k＜r)，每類 隨機(jī)抽取一個(gè)樣本Ci=?Si作為訓(xùn)練集中的樣本，訓(xùn)練集為從剩下的(r-k)個(gè)樣本中隨機(jī)抽取t個(gè)樣本作為測(cè)試集

第三步，模型建立。

訓(xùn)練集的k個(gè)樣本，其壓力系數(shù)為CT(n1×k)，幾何外形（縱坐標(biāo)）為yT(n2×k)，對(duì)壓力系數(shù)和幾何外形分別建立POD模型，命名為POD 1和POD 2，得到對(duì)應(yīng)的基模態(tài)系數(shù)矩陣YC(m1×r)和Yy(m2×r)，其中m1和m2滿 足m1?n1、m2?n2，其 大 小 通 過 分 析POD重構(gòu)誤差確定。此后，將YC和Yy分為作為BPNN的輸入和輸出來訓(xùn)練BPNN模型。其中，BPNN為三層，輸入層為m1維，輸出層為m2維，隱藏層神經(jīng)元個(gè)數(shù)為max{m1，m2}。

第四步，翼型反設(shè)計(jì)。

任意選取測(cè)試集中的樣本，其對(duì)應(yīng)的壓力系數(shù)為Cs。首先，通過POD 1模型得到壓力系數(shù)的基模態(tài)系數(shù)YC s；然后，帶入訓(xùn)練好的BPNN，得到預(yù)測(cè)幾何外形的基模態(tài)系數(shù)Yy s；最后，通過POD 2重構(gòu)得到預(yù)測(cè)的幾何外形ys。從而快速實(shí)現(xiàn)翼型反設(shè)計(jì)。

2 算例與分析

2.1 翼型庫(kù)建立

采用Hicks-Henne參數(shù)化方法，計(jì)算工況為亞聲 速 狀 態(tài)：Ma=0.35、α=2.79°、Re=6.5×106；跨聲速狀態(tài)：Ma=0.75、α=2.79°、Re=5.7×106，流 場(chǎng) 求 解 器 分 別 為Xfoil［13］和Fluent。以NA?CA0012翼型求解為例，對(duì)比Xfoil與Fluent得到的翼面壓力系數(shù)分布，如圖4所示，可以看出：兩者基本吻合，說明Xfoil在小攻角亞聲速計(jì)算狀態(tài)下具有較高的準(zhǔn)確性。而后，利用拉丁超立方試驗(yàn)方法生成500個(gè)翼型樣本，并通過流場(chǎng)求解器計(jì)算500個(gè)新翼型的壓力系數(shù)，得到壓力系數(shù)矩陣C。

圖4 NACA0012翼型翼面壓力系數(shù)分布對(duì)比Fig.4 Comparison of pressure coefficient distribution of NACA0012

2.2 聚類取樣

定義重構(gòu)（或預(yù)測(cè)）誤差為

式中：e1為均方根誤差；e2為均方誤差；分別為第i個(gè)坐標(biāo)位置的預(yù)測(cè)數(shù)據(jù)和真實(shí)數(shù)據(jù)，該數(shù)據(jù)在POD 1中表示壓力系數(shù)，在POD 2中表示幾何外形縱坐標(biāo)；N為總坐標(biāo)位置數(shù)目。

對(duì)壓力系數(shù)C隨機(jī)抽取90%進(jìn)行POD建模，10%進(jìn)行POD重構(gòu)誤差分析。不同基模態(tài)個(gè)數(shù)取值下壓力系數(shù)的重構(gòu)誤差如圖5所示，可以看出：可以發(fā)現(xiàn)，壓力系數(shù)的信息主要集中在前11階：基模態(tài)個(gè)數(shù)小于11時(shí)，重構(gòu)誤差隨著基模態(tài)個(gè)數(shù)增加減小得非?？欤划?dāng)基模態(tài)個(gè)數(shù)繼續(xù)增大時(shí)，重構(gòu)誤差的下降程度變緩，并逐漸接近0。當(dāng)基模態(tài) 個(gè) 數(shù) 取20時(shí)，亞 聲 速e1=3.7×10-3，e2=1.39×10-5，跨 聲 速e1=3.0×10-3，e2=9.97×10-6，誤差很小，因此壓力系數(shù)的降階維數(shù)s1設(shè)為20。

圖5 壓力系數(shù)的重構(gòu)誤差Fig.5 Reconstruction error of pressure coefficient

采用K-means算 法［18］對(duì)500組壓力系 數(shù)的基模態(tài)系數(shù)（20階）進(jìn)行聚類分析，聚類數(shù)目設(shè)置為200。根據(jù)1.4節(jié)的第二步，可以得到200個(gè)訓(xùn)練樣本。再在剩余樣本中隨機(jī)抽取10個(gè)作為模型的測(cè)試樣本。

2.3 模型訓(xùn)練

使用與2.2節(jié)相同的方法確定幾何外形的基模態(tài)個(gè)數(shù)取值。隨機(jī)抽取訓(xùn)練樣本的90%，并對(duì)其幾何外形建立POD 2模型，然后對(duì)另外10%進(jìn)行重構(gòu)誤差分析，得到幾何外形的重構(gòu)誤差如圖6所示。

圖6 幾何外形的重構(gòu)誤差Fig.6 Reconstruction error of geometry

從圖6可以看出：當(dāng)基模態(tài)個(gè)數(shù)為12時(shí)，亞/跨聲速的均方根誤差e1基本接近于0，說明前12階基模態(tài)可以高精度地重構(gòu)出幾何外形。

對(duì)比圖5和圖6，可以看出：跨聲速的重構(gòu)誤差小于亞聲速，這是由于在Hicks-Henne參數(shù)化中跨聲速的參數(shù)范圍小于亞聲速所致。

此時(shí)，確定了BPNN模型的輸入和輸出維度分別為20、12，取隱藏層神經(jīng)元個(gè)數(shù)為20，訓(xùn)練并建立從壓力系數(shù)的基模態(tài)系數(shù)到幾何外形的基模態(tài)系數(shù)的網(wǎng)絡(luò)。由于訓(xùn)練樣本量小（200個(gè)），且BPNN模型僅為三層，模型訓(xùn)練耗時(shí)很短。在In?tel i5-8500 CPU、16G RAM的PC上，聚類取樣耗時(shí)約為0.5 s，POD 1和POD 2的建模耗時(shí)均不到0.2 s，BPNN訓(xùn)練耗時(shí)為5～8 s，整個(gè)模型的建模時(shí)間相對(duì)于流場(chǎng)求解時(shí)間是很短的。

2.4 結(jié)果分析

在某計(jì)算工況下，給定目標(biāo)壓力分布，獲得滿足該壓力分布的翼型的過程稱為翼型反設(shè)計(jì)［20］。在亞/跨聲速下，分別對(duì)于所選取的10個(gè)測(cè)試樣本，將其壓力系數(shù)作為模型輸入，使用兩層POD+BPNN模型，預(yù)測(cè)出對(duì)應(yīng)的翼型，并與目標(biāo)翼型的外形進(jìn)行對(duì)比，得到均方根誤差分布如圖7所示。

圖7 翼型反設(shè)計(jì)的均方根誤差Fig.7 RMSEs of inverse-designed airfoils

取亞/跨聲速中誤差最大的第9個(gè)和第1個(gè)測(cè)試樣本，對(duì)比其目標(biāo)翼型與預(yù)測(cè)翼型，如圖8所示。對(duì)比預(yù)測(cè)翼型的壓力系數(shù)與目標(biāo)壓力系數(shù)，如圖9所示。

圖8 目標(biāo)翼型和預(yù)測(cè)翼型的對(duì)比Fig.8 Comparison of real and predicted airfoils

圖9 目標(biāo)和預(yù)測(cè)壓力系數(shù)的對(duì)比Fig.9 Comparison of real and target pressure coefficient

從圖8～圖9可以看出：預(yù)測(cè)翼型與目標(biāo)翼型無明顯差異，且預(yù)測(cè)翼型的壓力系數(shù)與目標(biāo)壓力系數(shù)曲線吻合較好，說明在亞/跨聲速工況下，兩層POD+BPNN模型的精度能夠滿足翼型反設(shè)計(jì)的要求。

同時(shí)，以亞聲速為例，將預(yù)測(cè)翼型與目標(biāo)翼型的縱坐標(biāo)做差，得到上下翼面翼型外形的預(yù)測(cè)誤差Δ1；將預(yù)測(cè)翼型的壓力系數(shù)與目標(biāo)壓力系數(shù)做差，得到上下翼面壓力系數(shù)的預(yù)測(cè)誤差Δ2。兩個(gè)誤差如圖10所示，可以看出：翼型縱坐標(biāo)的預(yù)測(cè)誤差的量級(jí)為10-4，而翼型縱坐標(biāo)的量級(jí)為10-2；壓力系數(shù)的預(yù)測(cè)誤差的量級(jí)為10-2，而壓力系數(shù)的量級(jí)為100，說明使用該模型進(jìn)行翼型反設(shè)計(jì)具有較高的準(zhǔn)確性，且預(yù)測(cè)翼型的翼面壓力系數(shù)分布與目標(biāo)壓力系數(shù)分布相一致。

圖10 翼型外形和壓力系數(shù)的預(yù)測(cè)誤差Fig.10 Prediction errors of airfoils and pressure coefficients

上述模型是由聚類取樣得到的200個(gè)樣本訓(xùn)練所得，為了進(jìn)一步分析訓(xùn)練樣本量對(duì)模型預(yù)測(cè)精度的影響，選取聚類數(shù)目k∈[20，490]，并隨機(jī)取10個(gè)樣本作為測(cè)試集，得到亞/跨聲速下、不同k下10個(gè)測(cè)試樣本的預(yù)測(cè)翼型外形與真實(shí)翼型外形的均方根誤差均值e1隨k的變化規(guī)律如圖11所示，可以看出：聚類數(shù)目小于120時(shí)，隨著聚類數(shù)目的增加，均方根誤差下降很快，說明訓(xùn)練樣本量過少不能夠表征出全部樣本的特征；當(dāng)聚類數(shù)目達(dá)到180以上，均方根誤差處于較低水平，繼續(xù)增加聚類數(shù)目不會(huì)使均方根誤差產(chǎn)生明顯變化，說明聚類取樣得到的180個(gè)樣本基本可表征出全部樣本的特征。本文選取200個(gè)樣本建模是合理的。

圖11 聚類數(shù)目對(duì)模型預(yù)測(cè)誤差的影響分析Fig.11 Influence of the of cluster number on the model error

2.5 超參數(shù)分析

BPNN的超參數(shù)包括網(wǎng)絡(luò)層數(shù)、隱藏層的神經(jīng)元個(gè)數(shù)、激活函數(shù)類型、訓(xùn)練函數(shù)類型、學(xué)習(xí)率等，這些超參數(shù)會(huì)影響模型的訓(xùn)練效率、預(yù)測(cè)精度。本文針對(duì)跨聲速算例，以隱藏層神經(jīng)元個(gè)數(shù)和激活函數(shù)為例，分析超參數(shù)對(duì)模型訓(xùn)練的影響。

（1）隱藏層神經(jīng)元個(gè)數(shù)

通過設(shè)置不同的隱藏層神經(jīng)元個(gè)數(shù)，分析模型的預(yù)測(cè)誤差和訓(xùn)練耗時(shí)，如圖12所示。

圖12 隱藏層神經(jīng)元個(gè)數(shù)對(duì)模型訓(xùn)練的影響分析Fig.12 Influence of the number of hidden layer neurons on model training

從圖12可以看出：當(dāng)神經(jīng)元個(gè)數(shù)小于10時(shí)，預(yù)測(cè)誤差較大，說明在本問題中，隱藏層神經(jīng)元個(gè)數(shù)過少不能夠完整地學(xué)習(xí)到樣本的特征；當(dāng)神經(jīng)元個(gè)數(shù)取值在10～50之間，預(yù)測(cè)誤差基本一致，但是訓(xùn)練耗時(shí)會(huì)隨著神經(jīng)元個(gè)數(shù)的增多呈現(xiàn)上升趨勢(shì)。因此，隱藏層神經(jīng)元個(gè)數(shù)的取值應(yīng)適當(dāng)。

（2）激活函數(shù)

選取6種不同的激活函數(shù)組合，分析不同組合下模型的訓(xùn)練耗時(shí)和預(yù)測(cè)誤差，如表2所示，可以看出：使用“tansig+purelin”和“l(fā)ogsig+purelin”均能保證較高的訓(xùn)練效率和預(yù)測(cè)精度。

表2 不同激活函數(shù)組合對(duì)模型訓(xùn)練的影響分析Table 2 Influence of the different activation functions on model training

3 結(jié)論

（1）本文發(fā)展了一種基于兩層POD+BPNN的翼型反設(shè)計(jì)方法，可以高效地建立從高維翼面壓力系數(shù)到高維翼型外形的映射。

（2）通過對(duì)壓力系數(shù)的基模態(tài)系數(shù)聚類來確定樣本集，可以在保證模型精度的前提下有效地降低訓(xùn)練樣本量，提高建模效率。

（3）針對(duì)亞/跨聲速算例，使用200個(gè)樣本建立的兩層POD+BPNN模型，其預(yù)測(cè)翼型外形的均方根誤差為10-4量級(jí)，對(duì)應(yīng)的壓力系數(shù)與目標(biāo)壓力系數(shù)吻合很好，說明該模型的精度可以滿足翼型反設(shè)計(jì)的要求。