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

張 宇，王曉亮

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

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

圖1 飛艇結(jié)構(gòu)示意

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

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

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

1.1 數(shù)值方法

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

(1)

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

(2)

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

(3)

1.2 網(wǎng)格無關(guān)性驗(yàn)證

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

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

圖4 網(wǎng)格無關(guān)性檢驗(yàn)

表1 邊界條件類型

圖2 計(jì)算域幾何示意圖

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

(4)

式中，V為飛艇體積.

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

(5)

Fd=Fp+Ff,

(6)

Cdv=Cdpv+Cdfv.

(7)

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

1.3 算法驗(yàn)證

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

圖5 數(shù)值方法驗(yàn)證結(jié)果

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

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

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

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

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

飛艇外形可由6階多項(xiàng)式表達(dá)：

(8)

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

(9)

3 全局敏感度分析方法

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

3.1 Sobol全局敏感度分析方法

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

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

(10)

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

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

(11)

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

(12)

積分式(11)有

(13)

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

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

(14)

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

(15)

式(15)右側(cè)第2項(xiàng)被稱為總方差，寫為

(16)

偏方差定義如下:

(17)

偏方差和總方差的關(guān)系為

(18)

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

(19)

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

D=Di+Dei+Di,ei,

(20)

引入新的參量STi，即

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

(21)

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

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

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

(22)

有關(guān)量可由蒙特卡洛方法近似得到[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之間，對(duì)應(yīng)艇身長細(xì)比處于2～10之間，能涵蓋大多數(shù)常規(guī)艇身；最大半徑位置xd處于0.3～0.5之間，以避免艇身最大截面位置過于靠前或靠后，以此減少艇身外形出現(xiàn)畸變的概率；尾部張開角βt處于20°～40°之間，以保證艇身尾部光滑收縮；尾部偏移角αt，尾部厚度Δyt和尾部高度yt被設(shè)定為0，以確保表征的外形為旋成體.如圖7所示，不恰當(dāng)?shù)膮?shù)組合會(huì)導(dǎo)致病態(tài)外形出現(xiàn)，實(shí)線外形包含了負(fù)值，虛線外形有多個(gè)極值，因此這些病態(tài)外形要在抽樣時(shí)被剔除.基于表3中的參數(shù)范圍，得到所有非病態(tài)外形樣本構(gòu)成的范圍如圖8所示，可見所取的參數(shù)能很好地保證艇身廓線的多樣性.

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

圖8 艇身廓線構(gòu)成的樣本空間

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

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

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

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

圖10 病態(tài)和設(shè)計(jì)空間

(27)

4.2 阻力特性分析

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

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

圖12 Cdv和xd的關(guān)系

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

5 結(jié) 論

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

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

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

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