參數(shù)化艇身阻力特性的全局敏感度及設計空間

張 宇，王曉亮

(上海交通大學 航空航天學院，上海 200240)

平流層一般指海拔高度處于10～55 km的空間，由于平流層的垂直對流效應小，氣流穩(wěn)定，適合布置平流層飛行器完成諸如中繼通信、監(jiān)察和運輸?shù)热蝿眨哂袠O大的軍事及民用價值[1].平流層飛艇是一種能定點飛行、效費比高的平流層飛行器，其結構一般包括艇身、吊艙、尾翼和支撐骨架(如圖1所示).對平流層飛艇而言，是否具備長航時是一項重要性能指標[2-3].平流層飛艇飛行過程中受到的阻力由動力系統(tǒng)平衡，阻力系數(shù)與動力消耗近似為正比關系，加之飛艇表面積巨大，駐留風阻很高，因此即便阻力系數(shù)略微減小，也能極大的減小動力系統(tǒng)的能量消耗，從而提高執(zhí)行任務時間.通常，艇身阻力占到了全艇阻力的約2/3[4-5]，因此找到使艇身具有最小阻力系數(shù)的外形具有十分重要的意義[6-7].

圖1 飛艇結構示意

目前，飛艇減阻一般采用優(yōu)化方法進行. Lutz等[8]通過在艇身布置周向分布的點源/匯從而得到了壓力場和速度場，并利用邊界層模型得到了不同雷諾數(shù)下的最小阻力外形.Wang等[9]采用勢流-邊界層耦合方法和混合遺傳算法對平流層飛艇艇身外形進行了優(yōu)化. Geruti等[10]基于粒子群優(yōu)化算法(PSO)提出了適用于考慮附加質量的非常規(guī)布局飛艇的優(yōu)化框架. Kanikdale等[11]采用GNVR作為飛艇基礎外形，提出了飛艇外形的多變量約束方法，并用模擬退火算法對外形進行了優(yōu)化. Wang等[12]將氣動、結構、能源和質量作為優(yōu)化對象，通過綜合目標方程對飛艇外形進行了多學科優(yōu)化.一般地，飛艇艇身外形廓線可由八參數(shù)、四參數(shù)、GNVR、橢球形和“海豚”形描述[13-15].但這些方法能包含的外形范圍較小，表達過于數(shù)學化，效率較低. 翼型參數(shù)化方法較上述飛艇外形描述方法豐富很多，能表征更為多樣的外形空間.翼型參數(shù)化方法主要有Hick-Henne Bump Functions、B-splines、PARSEC和Class/Shape function Transformation (CST)等4種.其中，Hick-Henne Bump Functions方法通過疊加形狀函數(shù)并改變形狀函數(shù)因子αi對已有的基礎外形進行擾動，例如: Zhong等[16]以RAE2822為基礎翼型，通過引入10個形狀函數(shù)來表示新的翼型；Yang等[17]利用開源軟件SU2，引入5個形狀函數(shù)分別對NACA0012和RAE2822翼型進行了優(yōu)化；總的來說，Hick-Henne Bump Functions方法適用于對已有外形進行優(yōu)化，難以形成設計空間.B-Splines方法通過在翼型周圍改變控制點位置來更新翼型形狀，例如Pérez等[18]使用B-spline方法控制螺旋槳不同展向處葉素的外形，通過設定不同的控制點數(shù)和容錯系數(shù)得到不同的三維螺旋槳CAD模型；Martin等[19]應用擴展的B-spline方法對三維機翼形狀進行了優(yōu)化；該方法具有很強的凸包性和穩(wěn)定性，但上述控制點沒有具體的物理含義，不便于形成明確的設計參考. CST方法通過類函數(shù)C(x)和形狀函數(shù)S(x)的乘積定義外形曲線，例如Masdari等[20]應用CST方法描述翼型形狀，結合離散渦方法，以功率系數(shù)為目標，對Savonius渦輪進行了優(yōu)化. 雖然該方法的適應性較強，但表達過于數(shù)學化，同樣不利于形成明確的設計參考. PARSEC方法通過11個參數(shù)控制翼型形狀，Chen等[21]通過PARSEC方法描述翼型形狀，以功率系數(shù)為目標對垂直軸流風力機進行了優(yōu)化，使功率系數(shù)比NACA0015翼型的結果提高了13.26%；PARSEC方法每個參數(shù)均有清晰的物理含義，利于形成明確的設計空間. Sripawadkul等[22]系統(tǒng)地研究了上述4種翼型參數(shù)化方法的特性，并給出了每種方法在不同指標下的評價結果，其中的正交性指標體現(xiàn)了各參數(shù)化方法中參數(shù)與外形是否一一對應，該指標在形成設計空間的過程中十分重要，結果表明只有PARSEC方法與CST方法能完全滿足正交性.考慮到CST方法表達過于繁復，可選PARSEC方法來描述飛艇外形.

一般而言，優(yōu)化過程是面向所有參數(shù)的，這往往導致在次要參數(shù)上會消耗許多不必要的時間.因此提取對艇身阻力系數(shù)最敏感的參數(shù)，降低優(yōu)化設計維度，并由這些參數(shù)形成飛艇艇身外形設計空間，從而給予艇身初期設計一定的參考是十分必要的.基于此，本文以艇身阻力系數(shù)為目標，通過PARSEC參數(shù)化方法、CFD方法和Sobol全局敏感度分析方法得到對阻力系數(shù)最敏感的參數(shù)，進而形成了由這些參數(shù)構成的飛艇艇身外形設計空間.

1 數(shù)值方法及其驗證

1.1 數(shù)值方法

艇身阻力系數(shù)可由勢流理論、風洞實驗和計算流體力學(CFD)方法得到[23].風洞實驗的高耗時和高成本不適用于本文工作，針對飛艇發(fā)展的勢流理論通常只對細長體預測的較準確，且黏性作用往往使勢流理論的結果不準確.隨著計算機技術的發(fā)展，CFD成為解決氣動優(yōu)化設計和流固耦合問題的重要手段.低速不可壓縮流動的NS方程為

(1)

式中：ρ為流體密度;p為壓力;μ為動力黏度;u為速度矢量.對式(1)進行無量綱化可得

(2)

選取Spalart-Allmaras模型求解湍流流動，Spalart-Allmaras模型適合求解航空外流場問題，其計算開銷低，穩(wěn)定性好[24].求解器采取不可壓縮形式的壓力基求解器.壓力-速度耦合格式為“coupled”，空間梯度離散方法為“l(fā)east squares cell based”，壓力項采用二階格式離散，動量和修正湍流黏度采用二階迎風格式離散.Spalart-Allmaras模型通過求解關于修正湍流黏度的輸運方程獲得流場信息:

(3)

1.2 網(wǎng)格無關性驗證

因艇身外形一般為旋成體，故簡化為二維軸對稱模型.艇身長度為L，艇身與入口和出口的距離均為60L，遠場邊界與旋轉軸的距離也為60L，如圖2所示.模型經(jīng)過歸一化處理，模型長度L取為1.邊界條件設置見表1.流場的數(shù)值求解精度主要由y+和艇身沿軸向的網(wǎng)格種子數(shù)量決定.設置艇身表面的第1層網(wǎng)格高度為10-6L，計算得到的流線型旋成體模型輪廓線上的y+分布如圖3所示，橫坐標為歸一化的艇身軸向位置，可見本文設置能使y+滿足要求.如圖4所示，對比了在不同Rev下艇身軸向網(wǎng)格種子數(shù)量與體積阻力系數(shù)Cdv的關系，可見當艇身軸向網(wǎng)格種子數(shù)量取為600時，Cdv幾乎不再變化.后續(xù)數(shù)值分析均基于以上述網(wǎng)格劃分模式進行.

圖3 不同Rev下Model 4154表面的y+分布

圖4 網(wǎng)格無關性檢驗

表1 邊界條件類型

圖2 計算域幾何示意圖

文中以飛艇體積的立方根值作為特征長度，因此Rev定義如下：

(4)

式中，V為飛艇體積.

總阻力系數(shù)Cdv由壓差阻力系數(shù)Cdpv和黏性阻力系數(shù)Cdfv構成：

(5)

Fd=Fp+Ff,

(6)

Cdv=Cdpv+Cdfv.

(7)

式中：Fd為總阻力;Fp為壓差阻力;Ff為黏性阻力.

1.3 算法驗證

文獻[25]中對于流線型旋成體做了詳細的流動實驗分析記錄，在此選取6組模型作為驗證對象，模型材質為紅木，長度統(tǒng)一為2.74 m，通過在水洞中拖曳模型獲取阻力，模型中軸線距離水面深度為2.74 m，模型標簽名分別為4 154、4 156、4 158、4 162、4 164和4 175，對應長細比分別為4、6、8、7、7和5.流體介質的密度和動力黏度分別為998.2 kg/m3和0.001 003 kg/m·s，確保在數(shù)值計算中使CFD與實驗的雷諾數(shù)相等. 對比結果如圖5所示，計算和實驗結果的平均相對誤差為1.5%，最大相對誤差為3.8%，滿足工程誤差許可，表明本文數(shù)值方法是可靠的.

圖5 數(shù)值方法驗證結果

為支撐數(shù)值方法中雷諾數(shù)能表征流動情況這一結論，現(xiàn)從具體算例上進行結果分析.以4 154模型為對象，取流動雷諾數(shù)Rev為5×106，固定來流速度U∞為10 m/s，通過調整流體密度ρ與流體動力黏度μ使Rev保持不變.各種參數(shù)配合下4 154模型受到的阻力情況見表2，可見當雷諾數(shù)Rev一定時，模型的阻力系數(shù)可視為常數(shù).因此不論從控制方程還是計算方法來看，雷諾數(shù)均能表征流動情況.

表2 Rev=5×106時的阻力系數(shù)

2 飛艇外形參數(shù)化方法

PARSEC方法通過11個參數(shù)控制翼型形狀，每個參數(shù)均有清晰的物理含義.鑒于飛艇艇身的旋成體特性，只需選取上半部分的參數(shù)，于是可通過8個參數(shù)描述飛艇外形，如圖6所示.

圖6 基于“PARSEC”的8參數(shù)飛艇艇身外形

飛艇外形可由6階多項式表達：

(8)

式中的待定系數(shù)可通過求解下式獲得[21]，xte為弦長：

(9)

3 全局敏感度分析方法

敏感度體現(xiàn)了變量自身的改變對系統(tǒng)的影響程度.一般地，敏感度分析方法分為局部敏感度分析和全局敏感度分析.例如，龍卷風圖,基于分化的方法和篩查法屬于前者;基于回歸的方法,基于方差的方法,轉換不變方法和基于密度的方法屬于后者[26].局部敏感度能體現(xiàn)輸入空間一點附近的特性，全局敏感度還考慮了參數(shù)之間的相互影響，結果更合理可靠.

3.1 Sobol全局敏感度分析方法

假設物理模型能被方程f(x)描述，輸入x=(x1,x2,…,xn) 是n維空間的一點且xi(i=1,2,…，n)遵循[0,1]上的均勻分布.全體x構成一個n維立方體:

Cn={x|0≤xi≤1;i=1,2,…,n}.

(10)

輸出項f(x)能被分解為

f1,2,…,n(x1,x2,…,xn),

(11)

式中f0為常數(shù)，且f1,2,…,s(x1,x2,…,xs)關于它自身變量的積分為0.

(12)

積分式(11)有

(13)

由于式(11)右側至少有一個變量是不同的，故具有正交性，再由式(12)可知:

0,(i1,i2,…,is)≠(j1,j2,…,jl).

(14)

對式(11)積分和平方,即

(15)

式(15)右側第2項被稱為總方差，寫為

(16)

偏方差定義如下:

(17)

偏方差和總方差的關系為

(18)

引入比率S1,2,…,s,即

(19)

對于具有n個輸入變量的離散樣本，本文需要計算每個S1,2…,s的值，但對于需要借助CFD軟件才能獲得的輸出，這樣勢必造成極大的時間開銷.基于此，本文將所有輸入變量分為兩個子集，子集Xi僅僅包含變量xi, 另一子集Xei包含除了xi的所有變量. 因此總方差可寫為

D=Di+Dei+Di,ei,

(20)

引入新的參量STi，即

STi=Si+Si,ei=1-Sei,

(21)

Si表征了變量xi單獨對總方差的影響，STi體現(xiàn)了變量xi對系統(tǒng)的“總體”影響，即不僅僅包含變量xi自身也包含了變量xi和余下變量的任意組合對輸出的影響.因此，本文稱Si和STi為變量xi的一階敏感度和全局敏感度.

3.2 離散數(shù)據(jù)的敏感度分析

事實上，本文經(jīng)常面對的是離散的輸入與輸出數(shù)據(jù)，其中沒有明確的函數(shù)關系.現(xiàn)在假設有一具有n個輸入變量的系統(tǒng)，抽樣產(chǎn)生N個樣本.計算之前，本文需要準備兩組輸入數(shù)據(jù)x1,x2,…,xn和x′1,x′2,…,x′n，寫成矩陣形式如下：

(22)

有關量可由蒙特卡洛方法近似得到[27]：

(23)

(24)

(25)

(26)

式中：fi,j為將變量xj替換為x′j所得到的輸出；f′i,j為將變量x′j替換為xj所得到的輸出.于是本文能通過上述算法和定義獲得輸入變量的一階或全局敏感度.

4 敏感度分析和艇身阻力特性研究

4.1 敏感度分析

表3給出了基于PARSEC方法的飛艇外形參數(shù)范圍.頭部半徑rh不超過艇身長度1；最大半徑rd處于0.05～0.25之間，對應艇身長細比處于2～10之間，能涵蓋大多數(shù)常規(guī)艇身；最大半徑位置xd處于0.3～0.5之間，以避免艇身最大截面位置過于靠前或靠后，以此減少艇身外形出現(xiàn)畸變的概率；尾部張開角βt處于20°～40°之間，以保證艇身尾部光滑收縮；尾部偏移角αt，尾部厚度Δyt和尾部高度yt被設定為0，以確保表征的外形為旋成體.如圖7所示，不恰當?shù)膮?shù)組合會導致病態(tài)外形出現(xiàn)，實線外形包含了負值，虛線外形有多個極值，因此這些病態(tài)外形要在抽樣時被剔除.基于表3中的參數(shù)范圍，得到所有非病態(tài)外形樣本構成的范圍如圖8所示，可見所取的參數(shù)能很好地保證艇身廓線的多樣性.

圖7 不合理的參數(shù)組合造成的病態(tài)外形

圖8 艇身廓線構成的樣本空間

表3 飛艇外形參數(shù)范圍

飛艇飛行高度設定為20 km，當?shù)乜諝饷芏葹?.088 9 kg/m3，空氣動力黏度為1.421 61×10-5kg/m·s，大氣壓強為5 529.31 Pa.一般地，平流層飛艇的繞流雷諾數(shù)在6.0×106和1.2×107之間，對一典型體積為3.0×105m3的飛艇而言，當來流速度為20 m/s時，其雷諾數(shù)為8.37×106.故本文取雷諾數(shù)為8.0×106進行后續(xù)研究，進行數(shù)值分析時調整流體密度使所有工況的繞流雷諾數(shù)一致.

圖9顯示了每個外形參數(shù)對各阻力系數(shù)的總敏感度.可以發(fā)現(xiàn)，參數(shù)rh和rd對總阻力系數(shù)的敏感度分別達到68.8%和66.0%，除了xd以外的其余參數(shù)對總阻力系數(shù)幾乎沒有影響，且xd的總敏感度僅為2.1%，故剩余參數(shù)可設為確定值.觀察由參數(shù)rh和rd組成的設計空間，該空間可近似看作被3條線段和1條三次曲線包圍形成，如圖10所示.該空間的數(shù)學表述為

圖9 不同參數(shù)的總敏感度

圖10 病態(tài)和設計空間

(27)

4.2 阻力特性分析

分別設置參數(shù)xd，y″d和βt的值為0.35，-0.5和30°.然后將設計空間離散，由同樣的計算條件得到參數(shù)rh和rd與各阻力系數(shù)的關系，如圖11所示，隨著設計空間向右上方推進壓差阻力系數(shù)逐漸升高，而黏性阻力系數(shù)與此趨勢恰好相反，二者的綜合效果使總阻力系數(shù)在設計空間上存在極小區(qū)域.由此可得到參數(shù)rh和rd的最佳組合，另外由圖12可知，在其他參數(shù)確定的情況下隨著xd的增加，總阻力系數(shù)先減小再增大，當xd=0.435時，總阻力系數(shù)取最小值.

圖11 rh和rd形成的關于各個阻力系數(shù)的云圖

圖12 Cdv和xd的關系

綜上所述，可獲得使艇身具有最小阻力系數(shù)的參數(shù)組合，該艇身外形的頭部半徑rh為0.012 2，最大半徑rd為0.95，最大半徑位置xd為0.435.上述所得艇身長細比為5.263，與文獻[28]中的結論吻合.在實際設計階段，頭部半徑rh的范圍可取為[0,0.06]，最大半徑rd的范圍可取為[0.08,0.12]，該范圍內的阻力系數(shù)大小浮動程度較小，均可較好的降低艇身阻力系數(shù).

5 結 論

1)對平流層飛艇而言，動力消耗近似正比于阻力系數(shù)，艇身阻力占到了全艇阻力的約2/3，在初期設計階段，快速找到最佳艇身外形具有十分重要的意義.為降低設計維度，有必要對外形參數(shù)進行全局敏感度分析，進一步所得設計空間可為類似設計提供參考.

2)對流線外形的旋成體，宜采取2維軸對稱模型和一方程Spalart-Allmaras湍流模型進行流場求解，所得阻力系數(shù)與實驗值吻合度較高，具有很高的精度.

3) 在8.0×106這一典型雷諾數(shù)下，平流層飛艇阻力對外形參數(shù)的敏感度由高到低依次為：rh，rd和xd，其他參數(shù)的敏感度相較于這3個量可忽略不計，因此評估飛艇氣動特性時應同時考慮以上3個參數(shù)，而不只考慮長細比這一參量.

4)從減阻的角度出發(fā)，在實際設計中頭部半徑rh的范圍可取為0～0.06，最大半徑rd的范圍可取為0.08～0.12.