李 多, 曹軍義, 張征宇, 黃敘輝(. 中國空氣動(dòng)力研究與發(fā)展中心 高速空氣動(dòng)力研究所, 四川 綿陽 6000；. 西安交通大學(xué) 機(jī)械學(xué)院, 陜西 西安 70049)

0 引 言

目前，風(fēng)洞試驗(yàn)主要采取傳統(tǒng)的單變量方法(One Factor at a Time, OFAT)開展，其原理是在一次試驗(yàn)中選擇一個(gè)自變量(如迎角)進(jìn)行改變，并保持其它自變量不變，通過測量響應(yīng)變量(如氣動(dòng)力和力矩)，來得到自變量和響應(yīng)變量的關(guān)系。但以美國為代表的航空航天發(fā)達(dá)國家發(fā)現(xiàn)，OFAT方法存在很多缺陷。例如一個(gè)自變量的改變必定會(huì)導(dǎo)致另外的自變量發(fā)生變化，如迎角在改變過程中會(huì)引起側(cè)滑角的改變，這是OFAT風(fēng)洞試驗(yàn)方法數(shù)據(jù)的誤差來源之一[1]。同時(shí)，OFAT方法還存在忽略了試驗(yàn)中各變量間的交互作用、試驗(yàn)結(jié)果表示不直觀以及需要采集的數(shù)據(jù)點(diǎn)較多而導(dǎo)致的試驗(yàn)成本高效率低等缺陷[2-4]。

1997年，美國航空航天局(NASA)蘭利研究中心開始倡導(dǎo)使用現(xiàn)代試驗(yàn)設(shè)計(jì)方法(Modern Design of Experiment，MDOE)來代替OFAT方法。MDOE方法的原理是：每次試驗(yàn)中改變多個(gè)變量的值，同時(shí)考察多個(gè)變量與氣動(dòng)力、力矩間的關(guān)系，其核心是響應(yīng)面法(Response Surface Method, RSM)。MDOE方法的優(yōu)勢是：1) 需要較少的數(shù)據(jù)點(diǎn)，成本低；2)提供了一個(gè)自變量影響響應(yīng)變量的更全面的觀察而不是普通數(shù)據(jù)圖的表示方法：響應(yīng)僅是一個(gè)自變量的函數(shù)，其余變量被固定在定值；3) 能給出除物理設(shè)置外的各自變量的響應(yīng)值；4) 不僅能精確估計(jì)出設(shè)計(jì)空間某個(gè)特定點(diǎn)上的響應(yīng)值，同時(shí)還能給出該響應(yīng)值的不確定度，從而顯著減小試驗(yàn)結(jié)果的不確定度；5) 能夠很好的表示出變量之間的交互作用；6) 不需要知道響應(yīng)曲面的具體數(shù)學(xué)形式[5]。

本文在我國2.4 m跨聲速風(fēng)洞中開展了MDOE試驗(yàn)，研究試驗(yàn)中涉及的試驗(yàn)區(qū)間確定、子區(qū)間劃分、數(shù)據(jù)點(diǎn)數(shù)量選取、數(shù)據(jù)點(diǎn)在子區(qū)間的分布等具體問題，驗(yàn)證MDOE方法在降低試驗(yàn)成本、優(yōu)化試驗(yàn)結(jié)果表達(dá)、體現(xiàn)變量間交互作用等方面的優(yōu)勢。

1 試驗(yàn)設(shè)計(jì)

選擇迎角a和側(cè)滑角β為兩個(gè)自變量，阻力系數(shù)CD為相應(yīng)的響應(yīng)變量，試驗(yàn)區(qū)間為迎角-4°～8°，側(cè)滑角-5°～5°。由于變量的個(gè)數(shù)為2，試驗(yàn)區(qū)間可以表示成一個(gè)笛卡爾坐標(biāo)系的平面，其中x軸代表迎角，y軸代表側(cè)滑角，區(qū)間中每個(gè)數(shù)據(jù)點(diǎn)表示迎角與側(cè)滑角的不同組合。

1.1 子區(qū)間劃分

為了避免響應(yīng)面模型中增加高階項(xiàng)，必須將試驗(yàn)區(qū)間劃分為不同的子區(qū)間，最后通過整合所有的毗鄰的變量子區(qū)間，就能夠得到一個(gè)精確且連續(xù)的響應(yīng)面方程[6]?？紤]到正負(fù)區(qū)間的不對稱因素及大迎角范圍非線性嚴(yán)重的特點(diǎn)，將迎角區(qū)間劃分為3個(gè)子區(qū)間，如圖1及表1所示。

圖1 子區(qū)間劃分Fig.1 Subspace divide

表1 子區(qū)間劃分Table 1 Subspace divide

1.2 數(shù)據(jù)量確定

數(shù)據(jù)點(diǎn)的個(gè)數(shù)決定了模型的精度與所采取方法的效率，必須精確計(jì)算。Box和Draper給出了一個(gè)計(jì)算每個(gè)子區(qū)間試驗(yàn)點(diǎn)數(shù)n的方程[7]:

(1)

式中，ta代表置信區(qū)間的半寬，查統(tǒng)計(jì)手冊為1.96，γ為試驗(yàn)精度，取0.0001，p代表回歸模型中參數(shù)的個(gè)數(shù)等于6，σ為數(shù)據(jù)標(biāo)準(zhǔn)差，以精度最高的滾轉(zhuǎn)力矩為例，其標(biāo)準(zhǔn)差不大于0.00007，代入式(1)中可求得:

(2)

因此在本試驗(yàn)中，每個(gè)子區(qū)間指定12個(gè)數(shù)據(jù)點(diǎn)，共需要36個(gè)數(shù)據(jù)點(diǎn)。

1.3 數(shù)據(jù)點(diǎn)分布

在各子區(qū)間內(nèi)，數(shù)據(jù)點(diǎn)的分布方式不同導(dǎo)致方差不同。本文采用中心復(fù)合試驗(yàn)設(shè)計(jì)方法(central composite design, CCD)分布各子區(qū)間的試驗(yàn)點(diǎn)[8]，其中獨(dú)立變量經(jīng)過線性編碼按-1、1來進(jìn)行分布。式(3)說明了轉(zhuǎn)換過程[9]。

(3)

圖2表示了編碼CCD方法中點(diǎn)的分布，除了四個(gè)拐點(diǎn)與星點(diǎn)外，重復(fù)性被指定在設(shè)計(jì)中心來抑制預(yù)測方差，由于每個(gè)子區(qū)間內(nèi)有12個(gè)數(shù)據(jù)點(diǎn)，所以在設(shè)計(jì)中心重復(fù)點(diǎn)的數(shù)量是4。

12個(gè)數(shù)據(jù)點(diǎn)為每個(gè)子區(qū)間提供12個(gè)自由度，其在子區(qū)間內(nèi)的分配方法如表2所示。

圖2 數(shù)據(jù)點(diǎn)分布Fig.2 Data points distribution

表2 自由度分配Table 2 Degree of freedom distribution

1.4 試驗(yàn)設(shè)計(jì)總結(jié)

圖3顯示了36個(gè)數(shù)據(jù)點(diǎn)在試驗(yàn)區(qū)間的分布，共有29個(gè)點(diǎn)，其中三個(gè)點(diǎn)重復(fù)了4次，兩個(gè)點(diǎn)重復(fù)了2次，其中重復(fù)2次的數(shù)據(jù)點(diǎn)用來確保各子區(qū)間響應(yīng)面模型的獨(dú)立性。

圖3 數(shù)據(jù)點(diǎn)在試驗(yàn)區(qū)間分布Fig.3 Data points distribution in experiment space

2 試驗(yàn)數(shù)據(jù)處理

2.1 響應(yīng)面模型建立

開展風(fēng)洞試驗(yàn)獲取36個(gè)數(shù)據(jù)點(diǎn)上的氣動(dòng)系數(shù)后，建立了阻力系數(shù)CD的響應(yīng)面模型，如圖4所示。

從圖4的響應(yīng)面模型中看出：采用MDOE試驗(yàn)方法后，每個(gè)子區(qū)間內(nèi)僅需12個(gè)數(shù)據(jù)點(diǎn)就能得到氣動(dòng)系數(shù)的響應(yīng)值，而OFAT方法則需要45個(gè)數(shù)據(jù)點(diǎn)，試驗(yàn)結(jié)束后耗氣量結(jié)果對比顯示，MDOE方法的耗氣量僅占OFAT方法的50%，降低了成本；MDOE風(fēng)洞試驗(yàn)方法提供了一個(gè)迎角和側(cè)滑角影響各氣動(dòng)系數(shù)的更全面的觀察而不是普通數(shù)據(jù)圖的表示方法：響應(yīng)僅是迎角的函數(shù)，側(cè)滑角被固定在定值；MDOE試驗(yàn)方法打破了單變量方法側(cè)滑角被物理固定在整數(shù)角度(2°、4°、6°等)的局限性，給出了整個(gè)設(shè)計(jì)區(qū)間的響應(yīng)值，在設(shè)計(jì)區(qū)間內(nèi)，迎角和側(cè)滑角可以以任意小數(shù)角度進(jìn)行組合；最后，MDOE試驗(yàn)方法很好的表示出了兩個(gè)變量迎角和側(cè)滑角之間的交互作用，以圖4(a)為例，若想減小阻力系數(shù)CD，在只減小迎角的情況下CD變化不明顯，但如果在減小迎角的同時(shí)增大側(cè)滑角，則CD的值會(huì)急劇下降。

氣動(dòng)系數(shù)在各子空間內(nèi)響應(yīng)面模型的多項(xiàng)式構(gòu)成如表3所示，多項(xiàng)式構(gòu)成中包含的項(xiàng)用“√”標(biāo)出。

(a) -4°～0°

(b) 0°～4°

(c) 4°～8°

表3 響應(yīng)面多項(xiàng)式構(gòu)成Table 3 Polynomial composition of response surface

2.2 方差分析

以0～4°迎角范圍內(nèi)阻力系數(shù)CD為例，說明響應(yīng)面回歸的過程，方差分析表如表4所示[10]。

首先觀察ANOVA表中的總效果，對應(yīng)回歸項(xiàng)的P值為0.000(即<0.001)，表明本模型總的來說是有效的，對應(yīng)失擬項(xiàng)的P值為0.512，其值遠(yuǎn)比臨界值0.05大，即表明本模型沒有失擬現(xiàn)象。

最后觀察各項(xiàng)效應(yīng)的顯著性，可以看出迎角、迎角*迎角、側(cè)滑角*側(cè)滑角三項(xiàng)高度顯著，而側(cè)滑角、迎角*側(cè)滑角兩項(xiàng)不顯著(P值>0.05)，應(yīng)將其刪除。

從改進(jìn)后的方差分析表5中可以看出，各項(xiàng)效應(yīng)均為顯著，同時(shí)與全模型方差分析表中的數(shù)據(jù)相比R-Sq(調(diào)整)的值略有上升而s的值略有下降，數(shù)據(jù)表明改進(jìn)后的模型比原模型更準(zhǔn)確。

表4 方差分析Table 4 Analysis of variance

表5 方差分析改進(jìn)Table 5 Secondary analysis of variance

2.3 模型驗(yàn)證

為了驗(yàn)證所建立響應(yīng)面模型的準(zhǔn)確性，在每個(gè)子區(qū)間β=0°的曲線上選取5個(gè)額外的試驗(yàn)點(diǎn)檢驗(yàn)其是否能夠落在擬合曲線的95%置信區(qū)間內(nèi)，檢驗(yàn)結(jié)果如圖5所示,其中藍(lán)色曲線代表上下置信區(qū)間，黑色曲線代表本試驗(yàn)擬合曲線，紅色點(diǎn)代表檢驗(yàn)點(diǎn)。

從圖5中看出，在每個(gè)子區(qū)間內(nèi)，中間范圍的點(diǎn)檢驗(yàn)結(jié)果較好，而邊界地區(qū)的點(diǎn)檢驗(yàn)結(jié)果較差，也反映了中心復(fù)合設(shè)計(jì)方法在中心區(qū)域誤差較小而邊界處誤差較大的特點(diǎn)。但基本上所有的點(diǎn)都能落入95%置信區(qū)間內(nèi)，說明響應(yīng)面模型精度足夠。

(a) -4°～0°

(b) 0°～4°

(c) 4°～8°

3 結(jié) 論

通過本次試驗(yàn)的開展，主要得到以下結(jié)論：

1) 對比得出，MDOE風(fēng)洞試驗(yàn)方法僅需傳統(tǒng)單變量試驗(yàn)方法20%的測量數(shù)據(jù)量即可完成試驗(yàn)，所需的吹風(fēng)時(shí)間約占傳統(tǒng)單變量試驗(yàn)方法的50%，同時(shí)在試驗(yàn)結(jié)果表達(dá)、更多響應(yīng)值的提供、不確定度估計(jì)、交互作用體現(xiàn)等方面較單變量方法更有優(yōu)勢；

2) 方差分析結(jié)果表明，迎角對氣動(dòng)系數(shù)的影響比側(cè)滑角更為顯著，正確反映了該類飛行器的氣動(dòng)特性變化規(guī)律；

3) MDOE風(fēng)洞試驗(yàn)方法能夠以氣動(dòng)系數(shù)響應(yīng)面模型和方差分析表的形式提供更為詳細(xì)、全面的風(fēng)洞試驗(yàn)數(shù)據(jù)，在飛行器風(fēng)洞試驗(yàn)的方案設(shè)計(jì)及數(shù)據(jù)分析中具有較為廣闊的應(yīng)用前景。

[1]張江, 秦永明, 馬漢東. 亞跨超風(fēng)洞現(xiàn)代試驗(yàn)設(shè)計(jì)方法研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2015, 33(3): 384-391

[2]戰(zhàn)培國. MDOE風(fēng)洞試驗(yàn)方法研究[J]. 航空科學(xué)技術(shù), 2011(6): 11-14

[3]張江, 秦永明, 馬漢東. 風(fēng)洞 MDOE 的形式實(shí)驗(yàn)設(shè)計(jì)方法研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2016, 34(1): 59-69

[4]唐志共, 王文正, 陳功, 等. 氣動(dòng)模型在現(xiàn)代氣動(dòng)試驗(yàn)設(shè)計(jì)中的應(yīng)用研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(2): 172-176

[5]Kaymaz I, McMahon C A. A response surface method based on weighted regression for structural reliability analysis[J]. Probabilistic Engineering Mechanics, 2005, 20(1): 11-17

[6]Garrent J H. Where and why article neural networks are applicable in civil engineering[J]. Journal of Computing in Civil Engineering, ASCE, 1994, 8(2): 129-130

[7]Erickson G E. Analysis of a split-plot experimental design applied to a low-speed wind tunnel investigation[R]. NASA/TM-20130218013

[8]De Loach R. Applications of modern experiment design to wind tunnel testing at NASA langley research center[C]//36thAIAA Aerospace Sciences Meeting&Exhibit, Reno, Nevada, January 1998, AIAA-1998-0713

[9]Landman D, Simpson J, Vicroy D, et al. Efficient methods for complex aircraft configuration aerodynamic characterization using response surface methodologies[C]//44thAIAA Aerospace Sciences Meeting & Exhibit, Reno, Nevada, January 2006, AIAA-2006-0922

[10]馬逢時(shí), 周暐. 六西格瑪管理統(tǒng)計(jì)指南: MINITAB使用指導(dǎo)[M]. 北京: 中國人民大學(xué)出版社, 2013.