基于數論網格法與Morris法的翼型氣動特性敏感性分析

文澤軍，孟祥恒，肖 釗，張 帆

（湖南科技大學機電工程學院，湖南 湘潭 411201）

翼型氣動特性是風力機設計的決定性因素［1］。通過改變翼型的厚度、彎度和前緣半徑等設計變量，可以設計出不同幾何外形和氣動性能的風力機翼型［2］。目前，CST（class function/shape function transformation，類型函數/形狀函數轉換）方法被廣泛應用于風力機翼型設計中，提高了翼型設計的效率［3］。在CST方法的應用過程中，參數化模型的設計變量通常在一定區間內變化，且變化帶有一定的隨機性，常被認為是服從正態分布的隨機變量［4］。該隨機變量會導致所設計的風力機翼型出現幾何偏差，從而可能導致風力機翼型的氣動性能不能滿足設計要求［5］。如Campobasso等［6］指出風力機翼型在設計及制造時隨機產生的幾何誤差會對風力機的功率特性造成一定影響。因此，在風力機翼型氣動設計中須考慮幾何不確定性的影響。國內外學者圍繞翼型氣動特性開展了大量研究。如：Ju等［7］基于蒙特卡洛模擬（Monte Carlo simulation，MCS）研究了風力機翼型氣動特性對翼型前緣的敏感性，但MCS法的樣本收斂性易受到模擬次數及樣本數的限制；Ernst等［8］采用拉丁超立方法研究了風力機翼型的相對厚度和相對彎度等幾何結構參數對氣動特性的影響，但拉丁超立方法存在樣本空間收斂性不足、不適用于高維運算等問題；王曉東等［9-10］研究了多項式混沌方法在隨機方腔流動模擬中的應用；姜昌偉等［11］采用嵌入式多項式混沌展開法研究了在隨機邊界下流動與傳熱問題的不確定性，但由于須對系統內部方程進行擴維處理，容易引入數值誤差；劉全等［12］研究了非嵌入式多項式混沌方法在拉氏計算中的應用，但在正交多項式基的構建中易出現矩匹配方程的奇異解，難以保證計算精度；為了較好地求解圓筒實驗模型，梁霄等［13］提出了一種多元非嵌入多項式混沌方法，但在預測非嵌入樣本點處的響應值時，樣本點的數量和位置都不是任意的，求解展開式系數時主觀性較大，容易引入近似誤差；鄔曉敬等［14］采用非嵌入多項式混沌方法開展翼型氣動特性全局敏感性分析，但因易受多項式模型的階數影響，近似解計算的主觀性較大；祝宏宇等［15］采用Sobol全局敏感性分析方法研究翼型氣動特性的敏感性，但該方法在求解高維運算問題時時間成本過高；Wang等［16］采用多點泰勒展開的近似伴隨方法分析翼型幾何的設計變量對氣動特性的影響，但近似伴隨方法依然易受到多點泰勒展開式項數的影響，導致樣本空間的計算精度不高；Peter等［17］對翼型氣動設計中的敏感性分析方法進行了總結。綜上可知，目前的研究方法存在計算效率低、樣本收斂性差及不適用于高維運算等不足，因此須對翼型氣動特性的敏感性作進一步研究。

近年來，數論網格（number theoretic net，NT-net）法被廣泛應用于工程統計分析領域。方開泰等［18］闡述了NT-net法在統計中的應用；劉飛等［19］闡述了NT-net法在極大似然估計中的應用，指出NT-net法具有較高的采樣效率和計算精度；賴斌等［20］闡述了NT-net法在流體隨機模擬中的應用，指出NT-net法可適用于高維隨機模擬。目前，將NT-net法應用于風力機翼型氣動特性敏感性分析的研究還不多見。

因此，本文以風力機翼型S832（以下簡稱翼型）為研究對象，基于CST方法構建翼型參數化模型，并提出一種基于NT-net法與Morris法的翼型氣動特性敏感性分析方法。首先，基于CST方法建立精度較高的翼型參數化模型；其次，應用NT-net法對參數化模型中的設計變量進行抽樣，并采用改進的Morris法分析NT-net法的抽樣結果，研究在翼型設計過程中各設計變量對翼型氣動特性的貢獻程度，并通過MCS進行驗證分析；最后，研究主要設計變量對翼型外形及氣動特性的影響，以期為翼型的設計提供理論參考。

1 翼型參數化模型及氣動特性

1.1 翼型參數化模型

1.1.1 翼型的幾何特征

翼型為風力機葉片沿展向的某一橫截面的形狀，其前緣呈圓弧形，后緣呈鈍形或尖形。弦線與前緣的交點為前緣點，弦線與后緣的交點為后緣點。其幾何特征與氣動特性密切相關。翼型的幾何特征如圖1所示。

圖1 翼型幾何特征示意Fig.1 Schematic of geometric characteristics of airfoil

1.1.2 翼型參數化方法

CST方法是根據翼型的幾何特征，通過類函數和型函數來表示翼型外形變化的參數化方法。其優勢是可以用少量的變量表達出可由計算機識別的翼型幾何形狀。因此，本文采用CST方法建立翼型參數化模型。以前緣點為原點O，以前緣點指向后緣點的方向為橫向即ψ向的正向，以垂直于弦線向上為縱向即ξ向的正向，建立坐標系O-ψξ，則翼面可以表示為：

式中：B(ψ )為類函數；G(ψ )為型函數；ΔLe為翼型后緣厚度；c為翼型的弦長。

其中：

式中：N1、N2為類函數系數，常取N1=0.5，N2=1.0；t為翼型Bernstein多項式階數，t=0，1，…，W；At+1為Bernstein多項式的待定系數，可由最小二乘法擬合得到。

參數化誤差Δy為：

式中：yψ為參數化翼型外形中橫坐標為ψ處對應的縱坐標。

1.2 翼型氣動特性

翼型氣動特性主要由升力系數Cl和阻力系數Cd表征。升阻比ε定義為升力系數與阻力系數的比，即：

Xfoil軟件是一款交互式獨立翼型的設計與分析軟件。與其他CFD（computational fluid dynamics，計算流體動力學）模擬軟件相比，Xfoil軟件具有收斂快、效率高和易操作等優點。基于線性變換的高階面元法，可利用Xfoil軟件計算Cl和Cd等。

2 翼型氣動特性敏感性分析方法

2.1 NT-net法

數論方法（number theory method，NTM）的實質是在b維單位立方體Mb上找到一個均勻散布的點集［21］，該點集通常稱為NT-net，偏差較小的點集稱為好格子點（good lattice point，GLP）集合。

設 P={gk(k=1，…，n)}為 Mb上的一個點集。對于r∈Mb，令U(r，P)為P中滿足gk≤r的點數，則P的偏差D（n，P）為：

式中：|0，r|為邊長為r的矩形，R(|0，r|)為矩形|0，r|的體積空間。

確定最優系數(h1，…，hb)是獲得GLP集合的關鍵環節。當維數b≤18時，應用生成矢量產生GLP集合。

令生成矢量 (n；h1，…，hb)滿足最大公約數(n，hi)=1，i=1，…，b，1≤hi

令 qki滿 足 1≤qki≤n，則 具 有 最 小 偏 差 的 集 合{gk=(gk1，… ，gkb)} 即為(n；h1，···，hb)的GLP集合。

2.2 基于NT-net法抽樣的敏感性分析方法

抽樣是指將b維立方體Mb的GLP集合的點映射到給定分布上的過程。敏感性分析常通過MCS來抽樣模擬各種輸入參數的可變性［22］。NT-net法被稱為偽MCS法，因其良好的樣本收斂性而被廣泛應用。Morris法是由Morris提出的全局敏感性分析方法，其敏感性指標為變量元效應的均值和標準差，可用來評價系統輸入量對系統輸出量的影響［23］。

為量化每個Bernstein多項式待定系數的不確定性影響，令待定系數為服從正態分布N(μ，σ)的設計變量 Z，Z=(z1，z2，…，zb)T，Z 的聯合分布函數為F(zi)。根據樣本數n及維數b，通過查找生成矢量或采用分圓域方法得到b維立方體Mb的NT-net及GLP集合，則GLP集合中的每個點gki表征b維設計變量Z的一個樣本值，由此可以建立翼型參數化模型的設計空間。

建立敏感性模型 Q=f(Z)，其中 zi～N(μ，σ )。對于通過NT-net法抽樣得到的gki，在第v次抽樣的變量元效應Ei，v及其均值 μi與標準差Si為：

針對元效應會在非單調函數中同時出現正負值的問題，Campolongo等［24］將均值μi改進為新均值：

敏感性分析的首要目標是計算并選擇均值較高或標準差較低的變量元效應。因此，定義變異系數ζ為新的敏感性指標，其為變量元效應的均值與標準差的比值，即：

2.3 基于NT-net法抽樣的敏感性分析流程

翼型氣動特性敏感性分析的流程如圖2所示，具體為：

圖2 翼型氣動特性敏感性分析流程Fig.2 Sensitivity analysis flow of aerodynamic characteristics of airfoil

1）確定設計變量Z、樣本數n和維數b，確定設計變量的分布特征；

2）根據維數b的大小選定適用方法進行NT-net抽樣并生成對應設計變量Z的GLP集合；

3）從設計變量Z的GLP集合中抽取樣本點，根據式（1）至式（5）構建翼型參數化模型的設計空間，利用Xfoil軟件求解翼型參數化模型的氣動特性；

4）基于改進的Morris法，分析設計變量的變化對氣動特性的影響，根據式（12）至（16）計算翼型氣動特性的變異系數ζ，并找出影響翼型氣動特性的主要設計變量。

3 算例分析

3.1 翼型參數化模型的構建及氣動特性敏感性分析

3.1.1 翼型參數化模型的構建

選取風力機專用翼型S832為研究對象。為了得到擬合精度較高的翼型參數化模型，取Bernstein多項式階數W=4，5，6，采用CST方法分別對翼型的上、下翼面進行擬合并構建參數化模型。將翼型的幾何外形坐標導入翼型設計軟件Profili，可以更直觀地表示翼型的幾何外形特征。在Profili中，縱坐標y/c表示翼型幾何外形的縱坐標y相對于翼型弦長c的位置關系，橫坐標x/c表示翼型幾何外形的橫坐標x相對于翼型弦長c的位置關系。翼型模型的參數化誤差如圖3所示。

圖3 翼型模型的參數化誤差Fig.3 Parameterized error of airfoil model

由圖3可知，Bernstein多項式階數對CST方法的擬合精度影響較大。適當提高Bernstein多項式階數可有效提高CST方法的擬合精度，但高階Bernstein多項式容易導致擬合過程呈病態化。當W≥5時，參數化誤差最大為9.7×10-4，已達到精度要求，因此本文選取W=5。采用5階Bernstein時，翼型上、下翼面的參數化模型均有6個Bernstein多項式待定系數。以上翼面為例，A1、A2主要控制翼型前緣半徑的大小，A3、A4主要控制翼型相對厚度和相對彎度的大小，A5、A6主要控制翼型后緣厚度的大小。通過CST方法對翼型幾何特征進行參數化擬合，構建出計算機可識別的參數化模型，如圖4所示。

圖4 翼型參數化模型Fig.4 Parameterized model of airfoil

3.1.2 翼型氣動特性敏感性分析

采用Xfoil軟件求解翼型的氣動特性。計算條件為穩態不可壓縮氣流，雷諾數Re=1.0×106，來流攻角θ=5°。

令各設計變量的均值為其標準差的10倍，各概率密度分布函數不完全相同，并采用NT-net法對全部設計變量分別進行抽樣，構建各設計變量樣本空間。總變量維數b=12，應用生成矢量產生GLP集合。為便于快速計算，根據前文所述Bernstein多項式待定系數的影響機制及翼型幾何外形的平順性要求，將變量維數分為4組，每組維數為3，查表得到生成矢量為（101；1，40，85），則GLP集合為{gk=(gk1，gk2，gk3) 。}

采用改進的Morris法進行翼型氣動特性的全局敏感性分析。以θ=5°時S832翼型的升阻比ε為輸出量，設計變量zi（i=1，2，…，12）為輸入量，其中z1，z2，…，z6和z7，z8，…，z12分別為依次對應的上、下翼面的設計變量。該12個設計變量是最小二乘法曲線擬合的數據點，依次控制5階Bernstein多項式中的12個待定系數。以升阻比對設計變量的變異系數ζ為敏感性指標，采用改進的Morris法進行全局敏感性分析。θ=5°時升阻比對各設計變量的變異系數如圖5所示。

圖5 θ=5°時升阻比對各設計變量的變異系數Fig.5 Variation coefficient of lift-drag ratio to design variables at θ=5°

由圖5可知，上、下翼面的升阻比對各設計變量的敏感性并不相同，其中對z9、z4、z1、z7、z6的敏感性較高。z4、z9控制待定系數A4、A9，間接控制翼型最大相對厚度和最大相對彎度的大小，表示翼型氣動特性主要受到翼型最大厚度和最大彎度的影響，其中以z9的變異系數最大，表明下翼面翼型最大相對厚度和最大相對彎度對翼型氣動特性的影響比上翼面翼型最大相對厚度和最大相對彎度對翼型氣動特性的影響大；z1、z7控制著待定系數A1、A7，間接控制翼型的前緣半徑大小，表示翼型氣動特性容易受到翼型前緣半徑的影響；z6控制待定系數A6，間接控制翼型后緣厚度大小，表示翼型后緣厚度也會影響翼型氣動特性，但比前緣半徑的影響程度要低。其余設計變量的變異系數均較小，表示升阻比對這些設計變量的敏感性較低，即這些設計變量對翼型氣動特性的影響程度較低。

3.2 MCS驗證

采用MCS法，每次抽取101個樣本，共計12次，并建立翼型參數化模型，調用Xfoil軟件求解翼型氣動特性，然后采用Morris法進行全局敏感性分析。相對于抽樣和建模時間，調用Xfoil軟件及進行全局敏感性分析的時間可以忽略。在同等條件下，采用MCS法和NT-net法進行1次抽樣的時間分別為205 s和180 s。采用MCS法進行全局敏感性分析所需的樣本計算總數為n（i+2）［25］，其中n=101，i=12，則采用MCS法進行全局敏感性分析的計算次數為101×（12+2）=1 414，而NT-net法的計算次數為 101×12=1 212。MCS法和NT-net法的建模時長及進行全局敏感性分析的計算次數如表1所示。

表1 MCS法和NT-net法的建模時長及進行全局敏感性分析的計算次數Table 1 Modeling time of MCS method and NT-net method and calculation times of global sensitivity analysis

由表1可知，相比于MCS法，采用NT-net法的建模時間減少了300 s，總計算次數減少了202次。因此，采用NT-net法可以有效提高翼型氣動特性敏感性分析的效率。

3.3 主要設計變量對翼型外形及氣動特性的影響

在來流攻角θ=5°的工況下，翼型的氣動特性主要受翼型幾何特征的影響。在基于NT-net法得到的各設計變量樣本空間中，分析主要設計變量對翼型幾何外形及氣動特性的影響。當設計變量變化時，翼型的最大相對厚度和最大相對彎度幾乎同時發生變化，因此本文不對翼型最大相對厚度和最大相對彎度進行單獨分析。后緣厚度對翼型氣動特性的影響比前緣半徑的影響小。因此，本文僅對前緣半徑、最大相對厚度和最大相對彎度對翼型外形及氣動特性的影響作深入分析。

3.3.1 前緣半徑對翼型外形及氣動特性的影響

選取樣本空間中間接控制上、下翼面前緣半徑大小的設計變量的最大值z1max、最小值z7min得到對應的待定系數A1、A7，其余待定系數與翼型S832一致，擬合出光滑的幾何外形，命名為翼型1，如圖6所示。

圖6 以z1max、z7min間接擬合的翼型1的幾何外形Fig.6 Geometric profile of airfoil 1 indirectly fitted with z1max、z7min

由圖 6可知：當z1、z7分別達到z1max、z7min時，翼型上、下翼面對應數據點間的距離比S832的大，所控制的待定系數的變化幅度也比翼型S832的大，因此翼型1的前緣半徑比翼型S832的大；針對前緣區域，上翼面的變化幅度比下翼面的變化幅度大；z1、z7對翼型其余部位的影響程度較小。

將翼型1的幾何外形坐標導入翼型設計軟件Profili和Xfoil，得到翼型1的前緣半徑和升阻比。其與翼型S832的對比如表2所示。

表2 翼型1與翼型S832前緣半徑和升阻比的對比Table 2 Comparison of leading edge radius and lift-drag ratio between airfoil 1 and airfoil S832

由表2可知，翼型1的前緣半徑比S832增加了0.05 mm，升阻比提高了0.659 1。因此，當前緣半徑適當增大時，翼型的升阻比提高，氣動特性增強。因此，設計時翼型的前緣半徑可適當取較大的值。

3.3.2 最大相對厚度和最大相對彎度對翼型外形及氣動特性的影響

選取樣本空間中控制上、下翼面最大相對厚度和最大相對彎度的設計變量最大值z4max、最小值z9min得到對應的待定系數A4、A9，其余待定系數與翼型S832的一致，擬合出光滑的幾何外形，命名為翼型2，結果如圖7所示。

圖7 以z4max、z9min間接擬合的翼型2的幾何外形Fig.7 Geometric profile of airfoil 2 indirectly fitted with z4max、z9min

由圖7可知：當z4、z9分別達到為z4max、z9min時，翼型2上、下翼面對應數據點間的距離比翼型S832的大，所控制的待定系數的變化幅度也比翼型S832的大，因此翼型2的最大相對厚度和最大相對彎度比翼型S832的大；針對翼型的中部區域，下翼面最大相對厚度和最大相對彎度的變化幅度比上翼面的大；z4、z9對翼型其余部位的影響較小。

同理，可得到翼型2的最大相對厚度、最大相對彎度和升阻比。其與翼型S832的對比如表3所示。

表3 翼型2與翼型S832最大相對厚度、最大相對彎度和升阻比的對比Table 3 Comparison of maximum relative thickness,maximum relative camber and lift-drag ratio between airfoil 2 and airfoil S832

由表3可知，翼型2的最大相對厚度和最大相對彎度分別比S832增加了1.42%、2.47%，升阻比減小了1.096 4。因此，當最大相對厚度和最大相對彎度增大時，翼型的升阻比減小，氣動特性減弱，這是由于來流攻角較小的緣故。因此，在來流攻角較小的工況下，翼型最大相對厚度和最大相對彎度可取較小值。

4 結 論

1）分析了Bernstein多項式階數對CST參數化翼型設計中擬合精度的影響。在參數擬合過程中，適當增大Bernstein多項式的階數可有效提高CST方法的擬合精度。當階數W=5時，采用CST方法可以達到擬合精度的要求。

2）基于NT-net法構建了翼型參數化模型。闡述了NT-net法的計算原理，基于NT-net法抽樣并構建了翼型參數化模型的設計空間。經MCS法驗證，采用NT-net法可以有效提高翼型參數化表達過程中不確定性傳播的擬合精度及效率。

3）通過NT-net法與Morris法的結合，確定了影響翼型氣動特性的主要幾何因素并提出改進策略。影響翼型氣動特性的主要幾何因素依次為最大相對厚度和最大相對彎度、前緣半徑和后緣厚度等，其中當來流攻角較小時，前緣半徑適當取較大值，最大相對厚度和最大相對彎度取較小值，可以有效提高該工況下翼型的氣動特性。