李啟良,曹冠寧,李 璇,楊志剛,鐘立元
(1.同濟(jì)大學(xué) 上海地面交通工具風(fēng)洞中心,上海 201804;2.同濟(jì)大學(xué) 上海市地面交通工具空氣動(dòng)力與熱環(huán)境模擬重點(diǎn)實(shí)驗(yàn)室,上海 201804)
汽車的氣動(dòng)阻力由內(nèi)部和外部阻力組成。無論是傳統(tǒng)內(nèi)燃機(jī)汽車還是新型電動(dòng)汽車都需要外部氣流進(jìn)入前艙將發(fā)動(dòng)機(jī)、空調(diào)等熱量帶到車外。在滿足散熱性能的同時(shí),不可避免地產(chǎn)生內(nèi)流阻力。在汽車開發(fā)時(shí),如不考慮內(nèi)流直接進(jìn)行車身氣動(dòng)優(yōu)化,可能得到優(yōu)化外形并不是真實(shí)狀態(tài)的最佳氣動(dòng)外形。因此有必要將內(nèi)流作為約束進(jìn)行車身氣動(dòng)優(yōu)化。
早期車身氣動(dòng)優(yōu)化主要基于道路試驗(yàn)和風(fēng)洞試驗(yàn)。隨著數(shù)值方法和計(jì)算機(jī)硬件的發(fā)展,數(shù)值仿真已全面應(yīng)用于車身氣動(dòng)優(yōu)化中。車身氣動(dòng)優(yōu)化包括單一參數(shù)的優(yōu)化和多參數(shù)優(yōu)化。單一參數(shù)優(yōu)化比較簡單和直觀。劉湘云等[1]通過數(shù)值仿真研究了前車窗角度變化對(duì)車身氣動(dòng)性能的影響,尋找到研究車型的最佳前車窗角度。車身是由多參數(shù)組成,各參數(shù)之間也會(huì)相互影響。為此研究人員從單一參數(shù)擴(kuò)展到兩個(gè)及以上參數(shù)。Khondge等[2]通過網(wǎng)格變形方法,對(duì)前、后風(fēng)窗傾角和側(cè)窗內(nèi)傾角進(jìn)行氣動(dòng)優(yōu)化,尋找到最佳傾角組合。Magazoni等[3]使用伴隨優(yōu)化方法進(jìn)行后視鏡外形優(yōu)化,獲得低阻后視鏡外形。郭建成等[4]針對(duì)轎車的前阻風(fēng)板的安裝位置和高度建立Kriging近似模型,通過多島遺傳算法進(jìn)行全局尋優(yōu),明顯改善了轎車的氣動(dòng)性能。容江磊等[5]將Kriging近似模型與帶精英策略的非支配排序遺傳算法相結(jié)合,對(duì)跑車尾翼的斷面形狀進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì),取得良好的優(yōu)化效果。
隨著油耗要求提高,迫切需要從整車角度進(jìn)行多參數(shù)優(yōu)化,以獲得更低阻力的車身外形。Han等[6]首先采用穩(wěn)態(tài)流動(dòng)計(jì)算和伴隨優(yōu)化方法對(duì)整車6個(gè)區(qū)域進(jìn)行氣動(dòng)優(yōu)化,獲得優(yōu)化外形,然后使用渦分離模型對(duì)該優(yōu)化外形進(jìn)行瞬態(tài)流場(chǎng)分析。Sun等[7]采用響應(yīng)面方法對(duì)早期開發(fā)階段的運(yùn)動(dòng)多功能汽車進(jìn)行氣動(dòng)外形優(yōu)化,得到比原始?xì)鈩?dòng)阻力系數(shù)低0.018的最優(yōu)外形。Lundberg等[8]通過神經(jīng)網(wǎng)絡(luò)和進(jìn)化優(yōu)化方法全自動(dòng)地進(jìn)行汽車外形優(yōu)化,得到最優(yōu)模型的氣動(dòng)阻力比原始模型降低13%。所在課題組在過去幾年開展了多參數(shù)建模和遺傳算法相結(jié)合的低阻車身氣動(dòng)優(yōu)化研究,得到氣動(dòng)阻力系數(shù)僅為0.129的低阻車身等一系列研究成果[9,10]。
回顧已有研究發(fā)現(xiàn),絕大多數(shù)的車身氣動(dòng)優(yōu)化都不含內(nèi)流,優(yōu)化外形可能并不是真實(shí)狀態(tài)的最佳外形。為此,本文結(jié)合了數(shù)值仿真和風(fēng)洞試驗(yàn),建立三廂轎車的內(nèi)流簡化模型和參數(shù)化模型,基于遺傳算法進(jìn)行氣動(dòng)優(yōu)化,獲得帶內(nèi)流約束的優(yōu)化外形。
選用手動(dòng)擋的三廂轎車作為研究對(duì)象,該車的長、寬和高分別為4521、1788和1492 mm。整車放置在氣動(dòng)-聲學(xué)整車風(fēng)洞的6分量氣動(dòng)天平中。通過天平測(cè)量整車受到的氣動(dòng)力,通過葉輪式風(fēng)速儀得到風(fēng)扇在不轉(zhuǎn)時(shí)通過散熱器的風(fēng)量,如圖1所示。試驗(yàn)速度從80 km/h增加到160 km/h,以20 km/h為間隔。所有試驗(yàn)均開啟邊界層抽吸和移動(dòng)帶,車輪處于旋轉(zhuǎn)狀態(tài)。
圖1 氣動(dòng)力和進(jìn)氣量測(cè)量Fig.1 Aerodynamic force and air input measurement
創(chuàng)建長、寬和高分別為9倍車長、10倍車寬和4倍車高的虛擬風(fēng)洞。整車模型放置在虛擬風(fēng)洞的地面上,車頭距離計(jì)算域的進(jìn)口3倍車長。
使用Hypermesh進(jìn)行計(jì)算域內(nèi)部件面網(wǎng)格的生成。車身外表面基于自身流動(dòng)特征合理布置面網(wǎng)格大小。后視鏡和A柱區(qū)域,網(wǎng)格大小為5 mm;在后風(fēng)窗和車尾的流動(dòng)分離區(qū)域,網(wǎng)格大小為6 mm;底盤和前后輪的網(wǎng)格大小分別為15 mm和10 mm;發(fā)動(dòng)機(jī)艙中發(fā)動(dòng)機(jī)、變速箱、冷卻部件等眾多零部件,網(wǎng)格大小為5 mm;其他區(qū)域,網(wǎng)格大小約為10 mm?;谝陨暇W(wǎng)格大小,整個(gè)計(jì)算域共創(chuàng)建319萬個(gè)面網(wǎng)格。
使用Star-CCM+軟件生成共2100萬個(gè)體網(wǎng)格。體網(wǎng)格類型為Trimmer網(wǎng)格,并在車身所有壁面和地面布置邊界層網(wǎng)格,以滿足非平衡壁面函數(shù)對(duì)y+的要求。為了能夠更精確地模擬整車流動(dòng)情況,在計(jì)算域內(nèi)分別設(shè)置了3個(gè)加密區(qū),圖2為計(jì)算域中間截面網(wǎng)格。
圖2 計(jì)算域中間截面網(wǎng)格Fig.2 Middle-section mesh of computational domain
為了建立三廂轎車參數(shù)化模型,有必要對(duì)內(nèi)流進(jìn)行簡化,從而更好地實(shí)現(xiàn)外形參數(shù)的優(yōu)化,減少單個(gè)樣本點(diǎn)的計(jì)算時(shí)間,提高優(yōu)化效率。內(nèi)流簡化需要保證簡化模型的進(jìn)氣量與實(shí)車相同。為此,建立了與真實(shí)車輛相同的開口面積,保留原有冷卻模塊,如圖3所示。從圖3可以看出,內(nèi)艙絕大多數(shù)部件被刪除,前端進(jìn)氣通過內(nèi)部流道經(jīng)冷卻模塊后直接流出到車底。
圖3 簡化模型Fig.3 Simplified model
仍使用Star-CCM+軟件對(duì)簡化內(nèi)流模型生成Trimmer體網(wǎng)格,體網(wǎng)格總數(shù)為1800萬個(gè)。體網(wǎng)格生成的所有設(shè)置與未簡化的整車模型相同。
計(jì)算域入口為速度入口,對(duì)應(yīng)車速為120 km/h。出口為壓力出口,設(shè)為0 Pa。兩個(gè)側(cè)面和頂面設(shè)為對(duì)稱面,地面設(shè)定為移動(dòng)壁面,車輪設(shè)置為旋轉(zhuǎn)壁面。冷卻模塊采用多孔介質(zhì)模型,其阻力特性由試驗(yàn)提供,冷凝器和散熱器的阻尼系數(shù),如表1所示。為了實(shí)現(xiàn)簡化模型的進(jìn)氣量與原來實(shí)車一致,需要改變冷卻模塊的阻尼系數(shù)來調(diào)整進(jìn)氣量。簡化模型氣動(dòng)計(jì)算方法與真實(shí)模型完全相同,只是冷卻模塊阻尼系數(shù)值不同。
表1 冷凝器和散熱器參數(shù)Table 1 Parameters of condenser and radiator
選用k-ε兩方程湍流模型[11],采用二階精度計(jì)算3000步后,殘差降到10-4并且監(jiān)測(cè)到氣動(dòng)阻力系數(shù)Cd值沒有明顯變化時(shí),可以認(rèn)為計(jì)算結(jié)果收斂。每個(gè)計(jì)算工況在24核的惠普刀片服務(wù)器上需計(jì)算24 h。本文約有20個(gè)計(jì)算工況,完成所有計(jì)算需480 h。
圖4 設(shè)計(jì)變量Fig.4 Design variable
基于建立的簡化內(nèi)流轎車模型,選取圖4所示的6個(gè)設(shè)計(jì)變量。結(jié)合原始參數(shù)和工程經(jīng)驗(yàn),確定如表2所示的參數(shù)變化范圍。
表2 優(yōu)化參數(shù)Table 2 Optimization parameters
利用CATIA軟件建模,根據(jù)6個(gè)設(shè)計(jì)變量確定車身上需要重建幾何模型的部分,在這些部分創(chuàng)建了628個(gè)控制點(diǎn)、127條樣條曲線、79個(gè)幾何曲面,使之與原始幾何形狀保持一致。利用樣條曲線捕捉描繪原始幾何模型的輪廓特征,由樣條曲線填充生成的曲面高度逼近原始幾何形狀,并允許曲面隨著樣條曲線形狀的改變而改變。樣條曲線的形狀由控制點(diǎn)控制,6個(gè)設(shè)計(jì)參數(shù)為獨(dú)立變量,控制點(diǎn)的坐標(biāo)是因變量,通過設(shè)計(jì)變量與因變量之間的映射最終建立了378個(gè)函數(shù)關(guān)系式。通過改變參數(shù)水平,使控制點(diǎn)的坐標(biāo)隨著函數(shù)關(guān)系式而相應(yīng)改變,從而改變樣條曲線的形狀,生成新的曲面造型,最終得到滿足不同參數(shù)水平的車身形狀,成為車身優(yōu)化的樣本點(diǎn)。圖5為重建后的幾何形狀。
圖5 重建后的幾何形狀Fig.5 Reconstructed geometry
對(duì)6個(gè)設(shè)計(jì)變量進(jìn)行參數(shù)組合,由于車型本身的限制,后擋風(fēng)玻璃傾角c與行李箱蓋傾角d的變化一致,只需要其中1個(gè)參數(shù)的變化值即可。對(duì)每個(gè)變量進(jìn)行選擇并組合,首先通過優(yōu)化拉丁方設(shè)計(jì)(Latin square design)[12]的方法確定參數(shù)組合,然后根據(jù)設(shè)計(jì)變量組合修改參數(shù)化模型得到多個(gè)模型。由于車身幾何較復(fù)雜,處理每個(gè)模型需要消耗大量的時(shí)間,受時(shí)間和計(jì)算資源的限制,第一代選取10組參數(shù)組合進(jìn)行計(jì)算。然后通過以下兩種遺傳算法對(duì)第一代的結(jié)果進(jìn)行優(yōu)化選擇。
(1)輪盤賭選擇法
(2)Taguchi方法
魯棒性參數(shù)設(shè)計(jì)(RPD)是一種獲得可控參數(shù)水平的方法,可用于將輸出平均值設(shè)置在期望目標(biāo)處,并使該目標(biāo)值的周圍變化達(dá)到最小值。
Taguchi方法提供了一些規(guī)則去簡化和標(biāo)準(zhǔn)化試驗(yàn)設(shè)計(jì),這個(gè)方法的關(guān)鍵是使用統(tǒng)計(jì)方法設(shè)計(jì)試驗(yàn),以找到整個(gè)群體的最佳參數(shù)。使用一組正交陣列設(shè)計(jì)試驗(yàn),并且同時(shí)進(jìn)行,此時(shí)使用正交陣列會(huì)顯著減少所需試驗(yàn)的數(shù)量。Taguchi方法使用信噪比(SNR)分析試驗(yàn)結(jié)果,SNR是信號(hào)變量與噪聲變量的比率。SNR分析的目的是確定變量的最佳組合以獲得最佳響應(yīng)量,其定義為:
SNR=-10logδ
(1)
由于每個(gè)參數(shù)水平下會(huì)有數(shù)個(gè)樣本,通過求取SNR,可以得到對(duì)于每個(gè)參數(shù)SNR隨參數(shù)水平的變化。根據(jù)SNR曲線,對(duì)于每個(gè)參數(shù)均可以挑選出的一個(gè)最佳水平。
在風(fēng)速為120 km/h時(shí),未簡化的整車模型的氣動(dòng)阻力系數(shù)的數(shù)值仿真值和試驗(yàn)值分別為0.358和0.335,誤差約為7%;它們的進(jìn)氣量分別為1.38 m3/s和1.30 m3/s,誤差約為6%。
經(jīng)多輪計(jì)算和調(diào)整,確定簡化模型冷卻模塊的阻尼系數(shù)為45.26 kg/m4,黏性阻尼系數(shù)為112.55×105kg/(m3s)。此時(shí)簡化模型的氣動(dòng)阻力系數(shù)為0.311,進(jìn)氣量為1.39 m3/s,實(shí)車模型的氣動(dòng)阻力系數(shù)為0.358,進(jìn)氣量為1.38 m3/s。之所以簡化模型的氣動(dòng)阻力系數(shù)比實(shí)車模型小,是因?yàn)榇蠖鄶?shù)內(nèi)艙部件被刪除,底盤拉平。
圖6 中截面速度云圖Fig.6 Mid-section velocity contour
圖6為簡化模型和實(shí)車模型的中截面速度云圖,兩圖中最大速度和最小速度的分布區(qū)域基本一樣。簡化模型的底盤部分被拉平,相較于實(shí)車模型,車身下方的速度會(huì)稍高。氣流流過車底,由于速度較高,使得簡化模型的車尾部的氣流滯留,低速區(qū)的范圍比實(shí)車模型稍大??傮w來看,簡化內(nèi)艙和底盤對(duì)車身周圍速度的影響較小。
從以上對(duì)比結(jié)果可見,雖然簡化后模型的氣動(dòng)阻力系數(shù)比實(shí)車模型小,但是簡化模型與實(shí)車模型的流動(dòng)結(jié)構(gòu)相差非常小,且簡化模型的進(jìn)氣量與實(shí)車模型一致,可使用簡化模型進(jìn)行多參數(shù)氣動(dòng)優(yōu)化。
表3給出了基于優(yōu)化拉丁方設(shè)計(jì)方法創(chuàng)建的第1代參數(shù)組合及計(jì)算結(jié)果。從第1代10組模型的計(jì)算結(jié)果中挑選出結(jié)果最好的5組參數(shù)進(jìn)行輪盤賭選擇,通過產(chǎn)生隨機(jī)數(shù),經(jīng)過25輪選擇后得到下一代5組參數(shù)組合并進(jìn)行計(jì)算,得到第2代的計(jì)算結(jié)果,見表4。由第2代5組模型的計(jì)算結(jié)果可以發(fā)現(xiàn):這5組模型優(yōu)化的效果并不明顯,氣動(dòng)阻力系數(shù)均下降得很緩慢,與上一代較好的模型相比,氣動(dòng)阻力并沒有明顯改善。
表3 第1代參數(shù)組合及計(jì)算結(jié)果Table 3 The first generation of parameter combinations and calculation results
表4 第2代參數(shù)組合及計(jì)算結(jié)果Table 4 The second generation of parameter combination and calculation results
利用Taguchi方法進(jìn)行優(yōu)化,基于第1代10組參數(shù)組合的仿真結(jié)果,由式(1)計(jì)算出每組的SNR值。Taguchi方法的結(jié)果由每個(gè)參數(shù)各個(gè)水平和SNR平均值的散點(diǎn)圖提供。每個(gè)參數(shù)水平的SNR平均值由相應(yīng)參數(shù)的每個(gè)水平的SNR之和除以該水平下相應(yīng)參數(shù)的數(shù)量而得到。根據(jù)圖7可以確定每個(gè)參數(shù)的最佳值,基于“響應(yīng)值越小、結(jié)果越理想”的原則,同時(shí)考慮方程為負(fù)相關(guān)函數(shù),確定最終每個(gè)參數(shù)的最大值為整車氣動(dòng)阻力系數(shù)的最佳水平。
SNR平均值的最大值對(duì)應(yīng)的設(shè)計(jì)參數(shù)組合即為最佳組合。根據(jù)該優(yōu)化方法得到的最佳參數(shù)組合重新修改幾何模型,并進(jìn)行仿真計(jì)算,表5給出了最佳組合和原始模型的參數(shù)組合及計(jì)算結(jié)果。由表5中可以看出,得到最佳組合的氣動(dòng)阻力系數(shù)為0.298,比原始值降低了4.2%。
通過以上遺傳算法的優(yōu)化過程,發(fā)現(xiàn)Taguchi方法優(yōu)化過程快速且效果也較好,故將Taguchi方法優(yōu)化的結(jié)果定為最佳優(yōu)化結(jié)果。圖8給出了最佳優(yōu)化模型和原始模型的各設(shè)計(jì)變量參數(shù)。由表5和圖8可以看出,優(yōu)化后的模型變化較明顯的是尾部行李箱蓋傾角和離去角,均有向內(nèi)收縮的趨勢(shì)。其他變量的傾角變化約為1°~2°。
下面分別利用Pareto圖(見圖9)和流場(chǎng)云圖(見圖10)對(duì)這兩個(gè)模型進(jìn)行分析。圖9為Pareto圖,目標(biāo)函數(shù)為氣動(dòng)阻力系數(shù),它受到各個(gè)設(shè)計(jì)變量的綜合影響。從圖9可以看出:車尾部離去角e對(duì)氣動(dòng)阻力系數(shù)的影響最大,并與該目標(biāo)函數(shù)呈現(xiàn)負(fù)相關(guān)關(guān)系。
圖7 各個(gè)參數(shù)與SNR的散點(diǎn)圖Fig.7 Scatter plot of each parameter and SNR
表5 原始模型與最佳模型的參數(shù)對(duì)比Table 5 Parameters comparison of the original model and the best model
圖8 各設(shè)計(jì)變量參數(shù)示意圖Fig.8 Schematic diagram of each design variable
圖9 目標(biāo)函數(shù)的Pareto分析Fig.9 Pareto analysis of objective function
圖10 尾部速度矢量圖Fig.10 Velocity vector of vehicle wake
除了單一變量對(duì)目標(biāo)函數(shù)有影響之外,還有“側(cè)風(fēng)窗傾角f-發(fā)動(dòng)機(jī)罩蓋傾角a”、“側(cè)風(fēng)窗傾角f-前擋風(fēng)玻璃傾角b”這些變量的交互作用對(duì)目標(biāo)函數(shù)的影響也較大,值得關(guān)注。
圖10為模型尾部的速度矢量圖,圖中用方框標(biāo)出了尾渦的位置。從圖中可以看出:從原始模型到優(yōu)化模型,上、下2個(gè)渦的位置均向后方移動(dòng)。渦流區(qū)中心的后移能夠適當(dāng)降低車尾受到尾渦的影響,提升了車尾壓力,降低了氣動(dòng)阻力系數(shù)。
(1)通過比較三廂轎車氣動(dòng)阻力系數(shù)和進(jìn)氣量的仿真與試驗(yàn)結(jié)果發(fā)現(xiàn),當(dāng)速度為120 km/h時(shí),兩者氣動(dòng)阻力系數(shù)分別為0.358和0.335,進(jìn)氣量分別為1.38和1.30 m3/s,兩者差異均較小,表明本文數(shù)值仿真方法是正確、可行的。
(2)建立了三廂轎車內(nèi)流簡化模型,它與實(shí)車模型前端進(jìn)氣量幾乎相等,流動(dòng)結(jié)構(gòu)基本一致,可用于建立多參數(shù)氣動(dòng)優(yōu)化模型,降低了建模難度,縮短了優(yōu)化時(shí)間,提升了優(yōu)化效率。
(3)建立了帶內(nèi)流的三廂轎車參數(shù)化模型,通過Taguchi方法得到最優(yōu)模型,其氣動(dòng)阻力系數(shù)為0.298。相比原始模型,氣動(dòng)阻力下降了4.2%。對(duì)比優(yōu)化模型與原始模型,發(fā)現(xiàn)車尾部的參數(shù)變化較大,處于車頭和車中部分的設(shè)計(jì)參數(shù)改變幅度較小。
參考文獻(xiàn):
[1] 劉湘云,馮俊虎,鄭智貞,等. 不同前車窗角度對(duì)汽車動(dòng)力性能影響的數(shù)值模擬[J]. 北京工業(yè)大學(xué)學(xué)報(bào),2009,35(7):987-990.
Liu Xiang-yun,Feng Jun-hu,Zheng Zhi-zhen,et al. Effect of different front window angles on the automobile dynamic performances[J]. Journal of Beijing University of Technology,2009,35(7):987-990.
[2] Khondge A,Sovani S,Verma G. Automation of vehicle aerodynamic shape exploration and optimization using integrated mesh morphing and CFD[C]∥SAE Paper,2011-01-0170.
[3] Magazoni F,Buscariolo F,Maruyama F,et al. Aerodynamic shape improvement for driver side view mirror of hatchback vehicle using adjoint optimization method[C]∥SAE Paper,2015-36-0156.
[4] 郭建成,谷正氣,容江磊,等. 基于Kriging模型的某轎車前輪阻風(fēng)板優(yōu)化[J]. 鄭州大學(xué)學(xué)報(bào):工學(xué)版,2011,32(3):125-128.
Guo Jian-cheng,Gu Zheng-qi,Rong Jiang-lei,et al. Optimization of air dam skirts of the front-wheels of a car based on kriging model[J]. Journal of Zhengzhou University(Engineering Science),2011,32(3):125-128.
[5] 容江磊,谷正氣,楊易,等. 基于Kriging模型的跑車尾翼斷面形狀的氣動(dòng)優(yōu)化[J]. 中國機(jī)械工程,2011,22(2):243-247.
Rong Jiang-lei,Gu Zheng-qi,Yang Yi,et al. Aerodynamic optimization of cross sectional shape for a sports car′s rear wing based on kriging surrogate model[J]. China Mechanical Engineering,2011,22(2):243-247.
[6] Han T,Kaushik S,Karbon K,et al. Adjoint-driven aerodynamic shape optimization based on a combination of steady state and transient flow solutions[C]∥SAE Paper,2016-01-1599.
[7] Sun S,Chang Y,Fu Q,et al. Aerodynamic shape optimization of an SUV in early development stage using a response surface method[C]∥SAE Paper,2014-01-2445.
[8] Lundberg A,Hamlin P,Shankar D,et al. Automated aerodynamic vehicle shape optimization using neural networks and evolutionary optimization[C]∥SAE Paper,2015-01-1548.
[9] 韋甘,楊志剛,李啟良. 低阻車身形體的參數(shù)化建模與氣動(dòng)試驗(yàn)[J]. 同濟(jì)大學(xué)學(xué)報(bào):自然科學(xué)版,2014,42(5):769-772,781.
Wei Gan,Yang Zhi-gang,Li Qi-liang. A parametric modeling method of low-drag car body and aerodynamic test[J]. Journal of Tongji University (Natural Science), 2014,42(5):769-772,781.
[10] 韋甘,楊志剛,李啟良. 基于分步遺傳算法的車身氣動(dòng)優(yōu)化[J]. 吉林大學(xué)學(xué)報(bào):工學(xué)版,2014,44(6):1578-1582.
Wei Gan,Yang Zhi-gang,Li Qi-liang. Aerodynamic optimization method for car body based on process costing genetic algorithm[J]. Journal of Jilin University(Engineering and Technology Edition),2014,44(6):1578-1582.
[11] Shih T H, Liou W W, Shabbir A, et al. A newk-εeddy viscosity model for high Reynolds number turbulent flows: model development and validation[J]. Computers and Fluids,1995,24(3):227-238.
[12] Mckay M D, Beckman R J, Conover W J. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics,1979,21(2):239-245.