基于NACA 翼型參數化方法的氣動優化設計

羅福星

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

0 引言

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

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

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

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

1 翼型流場分析

1.1 NACA 翼型參數化方法

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

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

圖1 翼型相關參數示意

中弧線表達式如下：

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

厚度表達式如下：

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

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

式中，（xu，yu）為上弧線坐標，（xd，yd）為下弧線坐標，θ為中弧線對應切線角度。

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

1.2 基于RANS 的數值模擬

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

RANS 所求解的連續方程和動量方程如下：

1.3 標準算例驗證

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

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

表1 求解格式選擇

表2 計算域網格設置

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

圖2 全計算域網格示意

圖3 翼型表面邊界層網格

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

表3 網格無關性驗證

圖4 翼型表面壓力系數對比

2 優化設計方法

采用基于代理模型的優化設計流程框架，以翼型的氣動特性參數作為目標函數開展優化分析，流程框架如圖5 所示。

圖5 翼型優化設計流程

2.1 響應面代理模型

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

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

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

2.2 NLPQL 優化算法

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

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

3 翼型的優化設計

3.1 二維翼型減阻優化

選擇NACA0012 作為優化的基準翼型，設計工況為：馬赫數取0.73，迎角取2.5°，雷諾數取6.5×106。

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

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

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

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

表4 設計變量限制范圍

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

圖6 三維響應面示意

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

3.2 優化結果對比

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

表5 優化前后參數對比

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

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

圖7 翼型優化外形對比

圖8 翼型優化壓力系數對比

對原始翼型和優化后翼型的流場計算結果進行基于Q 準則的渦識別，如圖9 所示，其中Q＞0 表示存在渦結構。可以看到設計工況下的流場主要存在以下兩方面區別：

圖9 基于Q 準則的渦識別對比

1）翼型在13%弦長位置增加向下的彎度后，上曲面前部下移，更貼近設計工況的氣流方向，使得優化后翼型上表面的渦區域更薄；

2）通過優化最大相對彎度及其位置，使翼型下表面后部渦尺度更小，且尾緣后方脫落渦位置延后，相較于原翼型，優化后翼型在尾緣較近區域內不存在明顯的渦結構。

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

表6 不同參數化方法結果對比

4 結論

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

（2）結合基于拉丁超立方的DOE、響應面代理模型和NLPQL 優化算法的優化設計框架可以開展較為高效的優化設計。

（3）與較為通用的基于CST 參數化方法的翼型優化方法相比，翼型優化設計框架可以在部分工況下得到效果相當的優化結果，體現了該框架具有一定的實用價值。