基于數(shù)論網(wǎng)格法與Morris法的翼型氣動特性敏感性分析

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

（湖南科技大學(xué)機(jī)電工程學(xué)院，湖南 湘潭 411201）

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

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

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

1 翼型參數(shù)化模型及氣動特性

1.1 翼型參數(shù)化模型

1.1.1 翼型的幾何特征

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

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

1.1.2 翼型參數(shù)化方法

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

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

其中：

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

參數(shù)化誤差Δy為：

式中：yψ為參數(shù)化翼型外形中橫坐標(biāo)為ψ處對應(yīng)的縱坐標(biāo)。

1.2 翼型氣動特性

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

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

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

2.1 NT-net法

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

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

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

確定最優(yōu)系數(shù)(h1，…，hb)是獲得GLP集合的關(guān)鍵環(huán)節(jié)。當(dāng)維數(shù)b≤18時，應(yīng)用生成矢量產(chǎn)生GLP集合。

令生成矢量 (n；h1，…，hb)滿足最大公約數(shù)(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來抽樣模擬各種輸入?yún)?shù)的可變性［22］。NT-net法被稱為偽MCS法，因其良好的樣本收斂性而被廣泛應(yīng)用。Morris法是由Morris提出的全局敏感性分析方法，其敏感性指標(biāo)為變量元效應(yīng)的均值和標(biāo)準(zhǔn)差，可用來評價系統(tǒng)輸入量對系統(tǒng)輸出量的影響［23］。

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

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

針對元效應(yīng)會在非單調(diào)函數(shù)中同時出現(xiàn)正負(fù)值的問題，Campolongo等［24］將均值μi改進(jìn)為新均值：

敏感性分析的首要目標(biāo)是計算并選擇均值較高或標(biāo)準(zhǔn)差較低的變量元效應(yīng)。因此，定義變異系數(shù)ζ為新的敏感性指標(biāo)，其為變量元效應(yīng)的均值與標(biāo)準(zhǔn)差的比值，即：

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

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

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

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

2）根據(jù)維數(shù)b的大小選定適用方法進(jìn)行NT-net抽樣并生成對應(yīng)設(shè)計變量Z的GLP集合；

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

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

3 算例分析

3.1 翼型參數(shù)化模型的構(gòu)建及氣動特性敏感性分析

3.1.1 翼型參數(shù)化模型的構(gòu)建

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

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

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

圖4 翼型參數(shù)化模型Fig.4 Parameterized model of airfoil

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

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

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

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

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

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

3.2 MCS驗證

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

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

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

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

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

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

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

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

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

將翼型1的幾何外形坐標(biāo)導(dǎo)入翼型設(shè)計軟件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。因此，當(dāng)前緣半徑適當(dāng)增大時，翼型的升阻比提高，氣動特性增強(qiáng)。因此，設(shè)計時翼型的前緣半徑可適當(dāng)取較大的值。

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

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

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

由圖7可知：當(dāng)z4、z9分別達(dá)到為z4max、z9min時，翼型2上、下翼面對應(yīng)數(shù)據(jù)點間的距離比翼型S832的大，所控制的待定系數(shù)的變化幅度也比翼型S832的大，因此翼型2的最大相對厚度和最大相對彎度比翼型S832的大；針對翼型的中部區(qū)域，下翼面最大相對厚度和最大相對彎度的變化幅度比上翼面的大；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。因此，當(dāng)最大相對厚度和最大相對彎度增大時，翼型的升阻比減小，氣動特性減弱，這是由于來流攻角較小的緣故。因此，在來流攻角較小的工況下，翼型最大相對厚度和最大相對彎度可取較小值。

4 結(jié) 論

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

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

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