基于NACA 翼型參數(shù)化方法的氣動優(yōu)化設(shè)計

羅福星

（國網(wǎng)西藏電力有限公司超高壓分公司，西藏 拉薩 850000）

0 引言

在飛行器機翼的設(shè)計過程中，高升阻比、低阻力等多方位的需求與約束，使得優(yōu)化設(shè)計需要對機翼氣動外形進行多變量、大范圍迭代。在此過程中，機翼幾何外形的參數(shù)化表達尤為重要，優(yōu)秀的參數(shù)化方法既可以滿足光滑的幾何外形，增加設(shè)計的魯棒性，又可以用較少的設(shè)計變量描述關(guān)鍵的翼型特性。常用的參數(shù)化方式有[1]：NACA 翼型、GA（W）翼型等由嚴(yán)格數(shù)學(xué)表達式給出翼型定義方式，由基準(zhǔn)翼型和擾動形函數(shù)線性疊加的形函數(shù)擾動法[2]，由Sobiesky[3]提出的特征參數(shù)描述法，利用正交基函數(shù)描述外形的正交基函數(shù)法[4]。此外，還有較為常用的CST 參數(shù)化方法，該方法由波音公司的Kulfan 等[5，6]提出，適應(yīng)性強。

而在氣動計算中，得益于發(fā)展成熟的算力和CFD 算法較好的魯棒性，RANS 計算方法在科研及工程過程中扮演著相當(dāng)重要的角色，相較于采用直接數(shù)值模擬（DNS）和大渦模擬（LES），其占用計算資源大幅減少。

采用較為成熟的NACA 四位數(shù)翼型參數(shù)化方式，翼型最大相對厚度、最大相對彎度、最大彎度位置的3 個參數(shù)作為設(shè)計變量，用以描述翼型外形。在氣動計算中采用RANS 計算方法，并結(jié)合SST-kw 湍流模型開展計算。

優(yōu)化流程中采用拉丁超立方法進行實驗設(shè)計（DOE），得到優(yōu)化設(shè)計的響應(yīng)面代理模型，在此基礎(chǔ)上依據(jù)二次拉格朗日非線性規(guī)劃（NLPQL）開展氣動優(yōu)化，由以上算法構(gòu)筑了一套機翼優(yōu)化設(shè)計方案。

1 翼型流場分析

1.1 NACA 翼型參數(shù)化方法

翼型參數(shù)化方法的選擇對氣動優(yōu)化設(shè)計中的結(jié)果優(yōu)劣及優(yōu)化效率有相當(dāng)關(guān)鍵的作用。采用典型的低速翼型——NACA 四位數(shù)翼型，該系列翼型為美國國家航空咨詢委員會提出，其應(yīng)用廣泛，在翼型發(fā)展史中有著重要地位。其為經(jīng)典的低速翼型，滿足低速飛機的較大升阻比、較低的最小阻力系數(shù)且低阻范圍寬以及失速過程緩和的要求。

NACA 四位數(shù)翼型的中弧線、厚度沿弦向的分布均由嚴(yán)格的數(shù)學(xué)表達式給出，翼型參數(shù)描述如圖1所示[7]。

圖1 翼型相關(guān)參數(shù)示意

中弧線表達式如下：

式中，h和hmax表示彎度及最大彎度，x為弦向位置，p為最大彎度位置，c為弦長。

厚度表達式如下：

式中t和tmax表示厚度和最大厚度。

上弧線和下弧線表達式為：

式中，（xu，yu）為上弧線坐標(biāo)，（xd，yd）為下弧線坐標(biāo)，θ為中弧線對應(yīng)切線角度。

將中弧線表達式（1）及厚度表達式（2）代入（3）式中即可得到NACA 四位數(shù)翼型的參數(shù)表達式。

1.2 基于RANS 的數(shù)值模擬

采用基于求解雷諾平均的N-S 方程（RANS）進行流場的數(shù)值模擬，其利用Reynolds 時均應(yīng)力簡化連續(xù)方程和動量方程。

RANS 所求解的連續(xù)方程和動量方程如下：

1.3 標(biāo)準(zhǔn)算例驗證

為驗證采用的氣動數(shù)值計算方法的準(zhǔn)確性，采用馬赫數(shù)為0.63、迎角為2 度工況下的NACA0012 翼型試驗數(shù)據(jù)[8]作為驗證依據(jù)。

流體域選擇流向往前15 倍、往后25 倍的特征長度作為計算域范圍，以減小遠場邊界的影響。采用RANS 求解，考慮計算工況選擇SST-kw 湍流模型，詳細計算格式見表1，網(wǎng)格及邊界條件設(shè)置見表2。

表1 求解格式選擇

表2 計算域網(wǎng)格設(shè)置

網(wǎng)格采用的是C-H 型拓撲構(gòu)建的結(jié)構(gòu)化網(wǎng)格，以使得網(wǎng)格離散方式更加貼合流場梯度的變化方向，且結(jié)構(gòu)化的網(wǎng)格節(jié)點編碼方式更有利于提高求解效率。網(wǎng)格具體如圖2～3 所示。

圖2 全計算域網(wǎng)格示意

圖3 翼型表面邊界層網(wǎng)格

用表2 的三套網(wǎng)格計算，結(jié)果見表3。表3 結(jié)果表明滿足網(wǎng)格無關(guān)性要求，為此后續(xù)計算選擇中等密度網(wǎng)格。獲取該算例的翼型表面壓力系數(shù)沿弦向的分布，如圖4 所示，可以看到CFD 的計算結(jié)果與試驗數(shù)據(jù)相差不大，表明采用的數(shù)值算法可靠。

表3 網(wǎng)格無關(guān)性驗證

圖4 翼型表面壓力系數(shù)對比

2 優(yōu)化設(shè)計方法

采用基于代理模型的優(yōu)化設(shè)計流程框架，以翼型的氣動特性參數(shù)作為目標(biāo)函數(shù)開展優(yōu)化分析，流程框架如圖5 所示。

圖5 翼型優(yōu)化設(shè)計流程

2.1 響應(yīng)面代理模型

在DOE 中通過計算足夠數(shù)量的采樣點構(gòu)建樣本空間，以此作為響應(yīng)面模型的輸入。樣本采用拉丁超立方設(shè)計，該方法原理為：在m維的設(shè)計空間中，將每一維的坐標(biāo)空間均分為n個子區(qū)間，在隨機取n個點時保證每個子區(qū)間只被研究一次，從而構(gòu)成樣本數(shù)為n的拉丁超立方采樣設(shè)計。在完成樣本空間的目標(biāo)函數(shù)計算后，可以據(jù)此構(gòu)建代理模型。

響應(yīng)面模型則是一種應(yīng)用廣泛、高魯棒性的代理模型，采用多項式函數(shù)來擬合設(shè)計空間，反應(yīng)設(shè)計變量和目標(biāo)函數(shù)的關(guān)系[9]。采用二項式多項式作為響應(yīng)面方程，其表達式如下：

此外，由于代理模型為近似模型，故需要取樣本空間中的部分設(shè)計點用于交叉驗證，進行代理模型的精度檢測。

2.2 NLPQL 優(yōu)化算法

在得到代理模型后，采用非線性序列二次規(guī)劃法（NLPQL）求解其最優(yōu)解。NLPQL 是應(yīng)用二級泰勒級數(shù)來展開目標(biāo)函數(shù)并將約束條件線性化，以此將非線性問題轉(zhuǎn)化為了二次規(guī)劃問題。該算法采用類牛頓矩陣Bk 定義Lagrange 函數(shù)的逼近和Hessian 矩陣的逼近，以此得到一個NLPQL 的子問題，表達式如下[10]：

式中：d為搜索方向，Bk為類牛頓矩陣，f（x）為目標(biāo)函數(shù)，g（x）為約束函數(shù)，q為變量，xu和xl分別表示邊界約束的上下限。

3 翼型的優(yōu)化設(shè)計

3.1 二維翼型減阻優(yōu)化

選擇NACA0012 作為優(yōu)化的基準(zhǔn)翼型，設(shè)計工況為：馬赫數(shù)取0.73，迎角取2.5°，雷諾數(shù)取6.5×106。

根據(jù)NACA 翼型參數(shù)化表達方式，以翼型最大相對厚度、最大相對彎度、最大彎度位置作為優(yōu)化的設(shè)計變量；而優(yōu)化設(shè)計的目標(biāo)函數(shù)為設(shè)計氣動工況下的阻力系數(shù)。此外，再考慮優(yōu)化問題的約束條件：

（1）幾何約束：最大相對厚度不小于原翼型；

（2）氣動約束：升力系數(shù)不低于原始翼型。

最終得到如（8）所示的優(yōu)化問題表達式。此外，補充設(shè)計變量的限制范圍，控制優(yōu)化翼型與原始翼型間的差異性，見表4。

表4 設(shè)計變量限制范圍

針對該減阻優(yōu)化問題，依據(jù)圖5 所示的優(yōu)化流程，首先采用拉丁超立方法進行DOE，構(gòu)建樣本空間，依次對設(shè)計點進行流體的數(shù)值仿真模擬，并依據(jù)約束條件篩選可行解集，依據(jù)該可行解集得到響應(yīng)面代理模型，其中關(guān)于最大彎度及最大彎度位置的三維響應(yīng)面如圖6 所示。

圖6 三維響應(yīng)面示意

在得到響應(yīng)面的基礎(chǔ)上，依據(jù)二次拉格朗日非線性規(guī)劃開展氣動優(yōu)化，該算例共計63 步迭代數(shù)，得到該代理模型下的最優(yōu)解，并對該設(shè)計點進行單獨的CFD 驗證，最終輸出優(yōu)化后翼型。

3.2 優(yōu)化結(jié)果對比

在該優(yōu)化框架下翼型優(yōu)化前后的參數(shù)見表5。從表5 可以看到，通過增加了向下的最大彎度，并將最大彎度位置后移，實現(xiàn)了在設(shè)計工況下翼型的阻力系數(shù)降低34%。

表5 優(yōu)化前后參數(shù)對比

此處最大厚度優(yōu)化前后無變化，是由于阻力基本與最大厚度正相關(guān)，而氣動約束驅(qū)使迭代主要朝著減小阻力的方向開展優(yōu)化，同時幾何約束限制了厚度的減小，故在該優(yōu)化算法下，最大厚度的初值即為優(yōu)化后的可行解。

優(yōu)化前后的外形及壓力系數(shù)對比見圖7 和圖8。由外形圖可以看到，因為NACA 翼型是根據(jù)中弧線疊加厚度的參數(shù)化表達方式，為此優(yōu)化后的翼型依舊保有較高的光順性，且優(yōu)化后翼型前緣位置的下翼面厚度大于上翼面厚度[11]。此外，優(yōu)化后的翼型具有較好的壓力系數(shù)分布規(guī)律，并且最大相對彎度、最大彎度位置的微小改變便能引起較大的阻力系數(shù)變化，說明阻力系數(shù)對外形變化相當(dāng)敏感。

圖7 翼型優(yōu)化外形對比

圖8 翼型優(yōu)化壓力系數(shù)對比

對原始翼型和優(yōu)化后翼型的流場計算結(jié)果進行基于Q 準(zhǔn)則的渦識別，如圖9 所示，其中Q＞0 表示存在渦結(jié)構(gòu)。可以看到設(shè)計工況下的流場主要存在以下兩方面區(qū)別：

圖9 基于Q 準(zhǔn)則的渦識別對比

1）翼型在13%弦長位置增加向下的彎度后，上曲面前部下移，更貼近設(shè)計工況的氣流方向，使得優(yōu)化后翼型上表面的渦區(qū)域更?。?/p>

2）通過優(yōu)化最大相對彎度及其位置，使翼型下表面后部渦尺度更小，且尾緣后方脫落渦位置延后，相較于原翼型，優(yōu)化后翼型在尾緣較近區(qū)域內(nèi)不存在明顯的渦結(jié)構(gòu)。

將優(yōu)化結(jié)果與參考文獻[12]中的結(jié)果進行對比見表6。從表6 結(jié)果可以看出，依據(jù)翼型參數(shù)化的表達方式能得到與CST 參數(shù)化方法效果相當(dāng)（在文獻中的對應(yīng)設(shè)計變量下）的優(yōu)化結(jié)果，而且僅采用3 個優(yōu)化設(shè)計變量，說明該參數(shù)化表達方式的變量能夠精準(zhǔn)地體現(xiàn)出翼型設(shè)計的關(guān)鍵變量，這與飛機設(shè)計手冊[13]中的表述規(guī)律一致，最大厚度在12%～17%時能得到較大的升力系數(shù)，而且最大彎度及最大彎度位置是影響阻力系數(shù)等氣動特性的關(guān)鍵參數(shù)。

表6 不同參數(shù)化方法結(jié)果對比

4 結(jié)論

（1）NACA 翼型參數(shù)化方法的設(shè)計變量少且變量有明顯的物理含義，參數(shù)的小幅變化能較大程度影響翼型的氣動特性，此外基于該表示函數(shù)的翼型光順性和連續(xù)性好，可以較為高效地構(gòu)建用于優(yōu)化設(shè)計的翼型，減少設(shè)計點數(shù)量，節(jié)省計算資源。

（2）結(jié)合基于拉丁超立方的DOE、響應(yīng)面代理模型和NLPQL 優(yōu)化算法的優(yōu)化設(shè)計框架可以開展較為高效的優(yōu)化設(shè)計。

（3）與較為通用的基于CST 參數(shù)化方法的翼型優(yōu)化方法相比，翼型優(yōu)化設(shè)計框架可以在部分工況下得到效果相當(dāng)?shù)膬?yōu)化結(jié)果，體現(xiàn)了該框架具有一定的實用價值。