固體燃?xì)鉁u輪火箭發(fā)動(dòng)機(jī)燃燒組織優(yōu)化研究①

李 江，王 偉，劉 洋，劉詩(shī)昌，楊 昀，楊 颯

(西北工業(yè)大學(xué)燃燒、熱結(jié)構(gòu)和內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

0 引言

固體燃?xì)鉁u輪沖壓發(fā)動(dòng)機(jī)(Solid Propellant Air Turbo Rocket，SP-ATR)有機(jī)結(jié)合了渦噴發(fā)動(dòng)機(jī)和固體沖壓發(fā)動(dòng)機(jī)2種單一推進(jìn)方式的工作特點(diǎn)，在臨近空間飛行器推進(jìn)領(lǐng)域得到較大發(fā)展，獲得了內(nèi)容廣泛的研究成果［1－4］。但也發(fā)現(xiàn)，驅(qū)渦推進(jìn)劑需對(duì)補(bǔ)燃室中摻混燃燒和做功能力進(jìn)行平衡，現(xiàn)有渦輪葉片材料性能又制約了高能推進(jìn)劑的使用，影響了SP-ATR性能的進(jìn)一步提升。

為克服以上問題，本課題組提出在推力需求較低時(shí)，直接選用潔凈碳?xì)渫七M(jìn)劑產(chǎn)生相對(duì)分子質(zhì)量低、固相顆粒少、溫度適當(dāng)?shù)母邏喝細(xì)怛?qū)動(dòng)渦輪做功;性能要求較高時(shí)，新增一個(gè)獨(dú)立含硼富燃燃?xì)獍l(fā)生器將高能燃?xì)庵苯右胙a(bǔ)燃室，參與空氣、驅(qū)渦燃?xì)獾膿交烊紵^程。

在以上2種工作方案中，驅(qū)渦燃?xì)馀c來流空氣的燃燒過程均為補(bǔ)燃室內(nèi)主要化學(xué)反應(yīng)過程，二者高效、均勻摻混，將直接決定發(fā)動(dòng)機(jī)性能優(yōu)劣。選用的碳?xì)渫七M(jìn)劑一次燃燒產(chǎn)物中的可燃成分以氣相為主，存在初始能量要求高、火焰?zhèn)鞑ニ俣鹊偷忍攸c(diǎn)，在氣態(tài)環(huán)境中不能單純依靠其與空氣間速度差形成的粘性剪切力提高穿透能力。因此，補(bǔ)燃室內(nèi)穩(wěn)定、高效二次燃燒組織技術(shù)成為影響發(fā)動(dòng)機(jī)性能的核心問題之一［5］。

在國(guó)外SP-ATR的研究過程中，僅對(duì)驅(qū)渦固體推進(jìn)劑的需求、選用標(biāo)準(zhǔn)、配方改進(jìn)等進(jìn)行論述，較少涉及固體碳?xì)渫七M(jìn)劑的二次燃燒技術(shù)。國(guó)內(nèi)的陳湘［6］、梁振欣［7］等的研究集中于整體性能的一維性能計(jì)算及以熱力計(jì)算為基礎(chǔ)進(jìn)行的參數(shù)影響規(guī)律研究，可借鑒依據(jù)相對(duì)較少。

本文在分析本課題組前期開展地面直連實(shí)驗(yàn)?zāi)承┕r補(bǔ)燃室未點(diǎn)燃問題，利用實(shí)驗(yàn)結(jié)果進(jìn)行精度校核基礎(chǔ)上，采用數(shù)值模擬方式進(jìn)行了SP-ATR典型工況適宜于燃燒組織方式初選，并由此結(jié)合替代模型優(yōu)化方法，進(jìn)一步開展了基于單級(jí)改進(jìn)V型支板的最優(yōu)構(gòu)型研究。

1 地面直連實(shí)驗(yàn)問題及分析

為研究SP-ATR摻混燃燒過程，隔離渦輪壓氣機(jī)組，利用多級(jí)喉道噴管進(jìn)行渦輪落壓過程模擬，設(shè)計(jì)獲得的直連實(shí)驗(yàn)發(fā)動(dòng)機(jī)構(gòu)型及結(jié)構(gòu)參數(shù)如圖1所示［8］。實(shí)驗(yàn)裝置包括驅(qū)渦燃?xì)獍l(fā)生器、進(jìn)氣裝置、輸運(yùn)管道、摻混段、補(bǔ)燃室、尾噴管等部件，也預(yù)留富燃燃?xì)獍l(fā)生器，以便開展其他相關(guān)實(shí)驗(yàn)。

圖1 SP-ATR直連實(shí)驗(yàn)發(fā)動(dòng)機(jī)Fig.1 Schematic diagram of direct-connect SP-ATR

現(xiàn)有研究已表明，為促進(jìn)驅(qū)渦碳?xì)渫七M(jìn)劑一次燃?xì)馀c來流空氣間穩(wěn)定、高效燃燒，需在補(bǔ)燃室中增加摻混燃燒裝置［9－10］。但SP-ATR適宜摻混器形式及其安裝方式仍不清晰，因此本課題組前期研究中，在摻混段、補(bǔ)燃室間分別增加圖2所示的六葉花、圓孔型摻混器，以強(qiáng)化驅(qū)渦燃?xì)?、來流空氣間的摻混燃燒過程。兩構(gòu)型在結(jié)構(gòu)上均由外部周向均布圓孔和位于突出圓弧面通氣結(jié)構(gòu)組成，在保證流通面積相同前提下，僅在通氣孔分布形式上存在差異。

高空工況模擬實(shí)驗(yàn)中部件工作正常，各參數(shù)達(dá)到設(shè)計(jì)要求，獲得了較真實(shí)可靠的流道參數(shù)數(shù)據(jù)，表明實(shí)驗(yàn)方法切實(shí)可行。但采用該方法進(jìn)行地面狀態(tài)模擬時(shí)，較優(yōu)當(dāng)量比條件下出現(xiàn)一次驅(qū)渦燃?xì)馀c空氣混合后直接排除，補(bǔ)燃室內(nèi)未發(fā)生二次燃燒現(xiàn)象。說明上述2種摻混器形式在燃燒組織較困難工況下強(qiáng)化摻混效果有限，仍需進(jìn)一步開展其他結(jié)構(gòu)形式燃燒組織方法研究與論證工作。

圖2 摻混器結(jié)構(gòu)圖Fig.2 Schematic of mixing-enhanced device

2 燃燒組織方案初選

考慮到數(shù)值模擬方式可較方便開展各類組織方式研究，且本課題組前期開展的地面直連實(shí)驗(yàn)研究已獲得部分可直接借鑒數(shù)據(jù)［8］，可用于指導(dǎo)計(jì)算方法參數(shù)確定。因此，本文采用該方式進(jìn)行方案初選。

2.1 數(shù)值方法

2.1.1 物理模型

綜合SP-ATR直連發(fā)動(dòng)機(jī)構(gòu)型特點(diǎn)、研究?jī)?nèi)容、計(jì)算效率要求，選用的計(jì)算模型僅由摻混段、摻混器、補(bǔ)燃室和尾噴管4部分構(gòu)成。其中，補(bǔ)燃室內(nèi)徑為φ130 mm、摻混段和補(bǔ)燃室段長(zhǎng)度為1 055 mm、尾噴管未達(dá)到完全膨脹，長(zhǎng)度63 mm?？紤]到摻混器結(jié)構(gòu)形式較復(fù)雜，計(jì)算區(qū)域內(nèi)采用非結(jié)構(gòu)網(wǎng)格劃分，并進(jìn)行針對(duì)性局部網(wǎng)格加密。

2.1.2 計(jì)算模型

采用有限體積法對(duì)下式所示的湍流、燃燒反應(yīng)的三維雷諾平均NavierSP-ATRStokes(NSP-ATRS)方程進(jìn)行求解［11］:

式中 W為守恒向量;向量F(W)、G(W)分別為無粘和粘性通量;向量H為體積力或化學(xué)反應(yīng)模型的源向量。

SP-ATR補(bǔ)燃室內(nèi)的實(shí)際流動(dòng)為三維非定常氣固兩相流動(dòng)過程，為簡(jiǎn)化計(jì)算，對(duì)流場(chǎng)作如下假設(shè):

(1)三維定常反應(yīng)流場(chǎng)，室內(nèi)氣體滿足理想氣體狀態(tài)方程p=ρRT;

(2)由實(shí)驗(yàn)中所選用碳?xì)渫七M(jìn)劑性質(zhì)、一次燃燒產(chǎn)物屬性、組分含量大小等因素，確定計(jì)算中燃?xì)饨M分僅含H2、CO、CH4三項(xiàng)氣相組成，對(duì)應(yīng)質(zhì)量分?jǐn)?shù)由熱力計(jì)算獲得，反應(yīng)為簡(jiǎn)單的一步總包反應(yīng);

(3)忽略重力等徹體力影響，補(bǔ)燃室內(nèi)流場(chǎng)與外界無熱量交換。

湍流模擬選用對(duì)自由剪切和分離流動(dòng)有較高精度的k-ω Menter SST模型［12］。氣相化學(xué)反應(yīng)采用渦耗散模型進(jìn)行，化學(xué)反應(yīng)如表1所示。

表1 氣相化學(xué)反應(yīng)方程Table 1 Reactions of gas phase

計(jì)算中涉及空氣進(jìn)口、燃?xì)馊肟凇⒐腆w壁面和出口邊界。其中，空氣、燃?xì)膺M(jìn)口均采用質(zhì)量流率邊界，給定質(zhì)量流率、溫度和各組分質(zhì)量分?jǐn)?shù);固體壁面為絕熱無滑移邊界;出口邊界為超音速壓強(qiáng)出口，邊界參數(shù)由二階外推獲得。

2.1.3 數(shù)值校驗(yàn)

考慮到SP-ATR二次燃燒數(shù)值計(jì)算相關(guān)參考依據(jù)的缺乏，文中選擇本課題組前期開展的表2所示的2種不同工況SP-ATR摻混燃燒實(shí)驗(yàn)結(jié)果，對(duì)所采用的數(shù)值方法進(jìn)行校驗(yàn)，以驗(yàn)證數(shù)值方法的精度和可靠性。表2中，m、T分別為質(zhì)量流量和溫度;下標(biāo)a、g分別表示來流空氣和驅(qū)渦燃?xì)狻?/p>

表2 實(shí)驗(yàn)工況Table 2 Operating conditions

2種工況下地面直連實(shí)驗(yàn)測(cè)得發(fā)動(dòng)機(jī)典型位置壁面壓強(qiáng)分布與采用文中數(shù)值方法獲得壓強(qiáng)曲線對(duì)比如圖3所示。可看出，2次實(shí)驗(yàn)中二者吻合程度都較理想，說明文中所確定的數(shù)值計(jì)算方法結(jié)果可靠，可較準(zhǔn)確地反映SP-ATR內(nèi)實(shí)際摻混、燃燒工作過程。

2.2 燃燒組織方式比較

SP-ATR摻混燃燒研究相對(duì)較少，缺乏可直接借鑒燃燒組織方式，本文參考固體沖壓發(fā)動(dòng)機(jī)、渦噴發(fā)動(dòng)機(jī)現(xiàn)有研究結(jié)論，確定燃燒組織方式初選方案，如表3所示。其中，進(jìn)氣方式、V型摻混器、波瓣摻混器結(jié)構(gòu)參數(shù)均采用現(xiàn)有研究推薦值［13－15］;環(huán)形摻混器結(jié)構(gòu)形式與圓孔、六葉花類似，其與波瓣摻混器結(jié)構(gòu)如圖4所示。

圖3 實(shí)驗(yàn)1、2測(cè)量值與計(jì)算壓強(qiáng)對(duì)比Fig.3 Comparison of experimental pressure distribution with computation in Case 1 and Case 2

表3 燃燒組織初選方案Table 3 Pilot study of mixing and combustion

圖4 環(huán)形和波瓣摻混器結(jié)構(gòu)圖Fig.4 Schematic of circular and lobed mixer

統(tǒng)一計(jì)算邊界條件為地面狀態(tài)模擬的實(shí)驗(yàn)1條件，采用上節(jié)確定數(shù)值方法，依次研究表3中各方案。為表征SP-ATR總體性能，定義發(fā)動(dòng)機(jī)比沖計(jì)算式為

式中 F、v分別為推力和速度;下標(biāo)e為噴管出口截面;其余參數(shù)與表2一致。

選擇相同工況下，直接摻混條件下采用文中數(shù)值方法獲得的比沖為基準(zhǔn)，得到圖5所示的各方案比沖變化?？煽闯?，在各自結(jié)構(gòu)參數(shù)推薦值條件下，6種方案對(duì)SP-ATR比沖影響存在較大差異。其中，單級(jí)V型和波瓣摻混器結(jié)構(gòu)較適宜。

與波瓣混流器相比，V型結(jié)構(gòu)簡(jiǎn)單，應(yīng)用范圍更廣，對(duì)SP-ATR燃燒組織更具借鑒意義。因此，本文針對(duì)其開展優(yōu)化研究，以進(jìn)一步提高發(fā)動(dòng)機(jī)性能。

圖5 方案初選結(jié)果Fig.5 Pilot study results of combustion mode

3 V型摻混器優(yōu)化研究

V型摻混器半角、混流器長(zhǎng)度、放置位置等參數(shù)耦合性較強(qiáng)，加之還需針對(duì)SP-ATR補(bǔ)燃室特點(diǎn)進(jìn)行適當(dāng)修改，加大了采用經(jīng)驗(yàn)參數(shù)尋優(yōu)難度和計(jì)算工作量。因此，文中結(jié)合數(shù)值方法基礎(chǔ)上，嘗試采用基于替代模型的漸進(jìn)優(yōu)化方法，以期快速準(zhǔn)確地獲得比沖性能較優(yōu)構(gòu)型方案。

3.1 替代模型構(gòu)造與優(yōu)化過程校驗(yàn)

Haupt函數(shù)為典型高度非線性、多極值點(diǎn)函數(shù)［16］，常用于復(fù)雜函數(shù)替代模型和優(yōu)化過程驗(yàn)證，表達(dá)式為

式中 x1∈［0，3.5］;x2∈［0，3.5］。

根據(jù)現(xiàn)有研究結(jié)論，本文采用精度較高的徑向基函數(shù)(Radial Basis Function，RBF)構(gòu)造近似模型［17］。圖6給出了利用試驗(yàn)設(shè)計(jì)方法［18］產(chǎn)生50樣本點(diǎn)構(gòu)建的近似模型與原函數(shù)計(jì)算值比較?？煽闯?，二者形狀吻合較好，替代模型較好體現(xiàn)了原函數(shù)中所有關(guān)鍵點(diǎn)特性。

圖6 Haupt函數(shù)計(jì)算與替代模型近似值比較Fig.6 Comparison between calculation and surrogate model for Haupt function

替代模型除要求較好體現(xiàn)原函數(shù)各局部特征外，還需驗(yàn)證其數(shù)值誤差。文中在定義域內(nèi)利用隨機(jī)函數(shù)分別產(chǎn)生10、100個(gè)驗(yàn)證點(diǎn)，對(duì)上述模型進(jìn)行精度分析，預(yù)測(cè)結(jié)果與實(shí)際值對(duì)比如圖7所示?？煽闯觯?次產(chǎn)生的驗(yàn)證點(diǎn)獲得的函數(shù)近似值與真實(shí)值間的差異均較小。經(jīng)進(jìn)一步統(tǒng)計(jì)，二者的總體誤差分別為 2.296%、5.148%。

圖7 Haupt替代模型誤差驗(yàn)證Fig.7 Errors validation of RBF model for Haupt function

采用在局部及全局極值點(diǎn)尋優(yōu)過程中具有較好效果的多島遺傳算法(Multi-Island Genetic Algorithm，MIGA)，計(jì)算Haupt函數(shù)極小值點(diǎn)，以完成替代模型優(yōu)化過程校驗(yàn)［19］。表4給出了利用50個(gè)樣本點(diǎn)構(gòu)建的原函數(shù)、替代模型，在定義域內(nèi)進(jìn)行1 000次迭代后獲得的優(yōu)化結(jié)果比較?？煽闯?，采用RBF模型獲得的函數(shù)極小值誤差約為5.362 6%，滿足精度要求。

表4 Haupt函數(shù)優(yōu)化驗(yàn)證Table 4 Optimization validation of RBF for Haupt function

綜上可認(rèn)為，文中所選擇替代模型優(yōu)化方法切實(shí)可行，結(jié)合上節(jié)數(shù)值模擬方法，可用于V型摻混器較優(yōu)構(gòu)型研究。

3.2 V型摻混器構(gòu)型優(yōu)化

根據(jù)V型摻混器已有研究結(jié)論及SP-ATR應(yīng)用特點(diǎn)，文中發(fā)展了圖8所示的燃燒增強(qiáng)裝置:頂角邊改由3段構(gòu)成，可增加調(diào)節(jié)靈活性;尾部增加一段直段，在影響軸線附近流動(dòng)的同時(shí)，可方便開展斜坡交錯(cuò)結(jié)構(gòu)、富燃燃?xì)饧尤敕绞降刃问降倪M(jìn)一步設(shè)計(jì)與驗(yàn)證。

耦合性較強(qiáng)參數(shù)數(shù)量較多時(shí)，不利于替代模型正確反映復(fù)雜模型內(nèi)部特征。因此，將V型摻混器放置于補(bǔ)燃室軸線處，選擇圖9所示的6個(gè)獨(dú)立變量:半角α1、α2、α3，長(zhǎng)度 lf、lb、摻混器頂點(diǎn)距燃?xì)馊肟谖恢?l為設(shè)計(jì)變量，其余參數(shù)在初步設(shè)計(jì)時(shí)取為定值。

圖8 V型摻混器Fig.8 Schematic of V-shape mixing-enhanced device

圖9 設(shè)計(jì)變量示意圖Fig.9 Design variables

由SP-ATR補(bǔ)燃室構(gòu)型特點(diǎn)，適當(dāng)擴(kuò)展各變量定義域，確定其取值范圍為 α1∈［10，20］，α2∈［20，30］，α3∈［30，50］，lf∈［20，100］，lb∈［20，100］，l∈［100，500］。利用試驗(yàn)設(shè)計(jì)法，在上述定義域內(nèi)產(chǎn)生49個(gè)設(shè)計(jì)點(diǎn)，以保證其代表性和均勻性，部分設(shè)計(jì)點(diǎn)如表5所示。采用上節(jié)數(shù)值方法，對(duì)以上各樣本點(diǎn)進(jìn)行逐一計(jì)算，構(gòu)建樣本數(shù)據(jù)庫(kù)。解耦后，得到圖10所示的各獨(dú)立變量與比沖間變化關(guān)系。

表5 部分設(shè)計(jì)點(diǎn)Table 5 Partial design points

圖10 各設(shè)計(jì)變量與比沖曲線Fig.10 Relationship between design variable and specific impulse

可看出，在文中所采用定義域內(nèi)，采用V型摻混器后，比沖性能都比直接摻混狀態(tài)有所提升。經(jīng)進(jìn)一步開展參數(shù)敏感性分析可發(fā)現(xiàn)，半角α2、α3、α1變化影響更明顯，其他3個(gè)參數(shù)敏感性較差。

以該工作狀態(tài)下SP-ATR比沖指標(biāo)最大為優(yōu)化目標(biāo)，定義其優(yōu)化模型為

表6給出了利用MIGA優(yōu)化算法在定義域內(nèi)進(jìn)行1 000次迭代尋優(yōu)、對(duì)應(yīng)條件下數(shù)值模擬獲得的3組結(jié)果對(duì)比。可看出，雖然在設(shè)計(jì)變量具體取值上存在一定差異，但比沖優(yōu)化結(jié)果較一致，且三者與對(duì)應(yīng)結(jié)構(gòu)下數(shù)值結(jié)果差異也較小，說明文中所用優(yōu)化方法可靠性和穩(wěn)定性較好。

表6 優(yōu)化結(jié)果Table 6 Optimization results

經(jīng)與直接摻混條件下SP-ATR性能對(duì)比，3組優(yōu)化值結(jié)構(gòu)下，比沖增幅約13.50%，表明采用改進(jìn)后的V型摻混器方案，可獲得性能較優(yōu)的SP-ATR。

4 結(jié)論

(1)文中采用的SP-ATR補(bǔ)燃室摻混燃燒數(shù)值方法，可較真實(shí)反映發(fā)動(dòng)機(jī)內(nèi)部真實(shí)流動(dòng)狀態(tài)，利用該法燃燒組織方式初選表明，單級(jí)V型摻混器有利于發(fā)動(dòng)機(jī)性能提高。

(2)基于替代模型漸進(jìn)優(yōu)化可靠性、穩(wěn)定性較理想，結(jié)合文中數(shù)值模擬方法，可較好實(shí)現(xiàn)SP-ATR摻混器優(yōu)化研究。

(3)典型工況下優(yōu)化結(jié)果顯示，采用改進(jìn)后V型摻混器的SP-ATR比沖性能較直接摻混條件提高約13.50%，性能較理想。

(4)地面狀態(tài)直連實(shí)驗(yàn)中，較優(yōu)當(dāng)量比條件下，圓孔、六葉花型摻混器強(qiáng)化摻混效果有限。

［1］Bossard J A，Christensen K L，Poth G E.ATR propulsion system design and vehicle integration［C］//24th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.Boston，Massachusetts.AIAA 88-3071.

［2］Thomas M E，Bossard J A，Ostrander M J.Addressing emerging tactical missile propulsion challenges with the solid propellant air-turbo-rocket［C］//36th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit.Huntsville，Alabama，2000.

［3］Christensen K L.Vehicle performance optimization utilizing the air turbo-ramjet propulsion system:methodology development and applications［D］.Rolla:University of Missouri，1997.

［4］史亞紅，戴耀松.國(guó)外戰(zhàn)術(shù)彈道空氣渦輪沖壓發(fā)動(dòng)機(jī)技術(shù)分析［C］//98年聯(lián)合推進(jìn)會(huì)議.北京，1998.

［5］李文龍，李平，郭海波.空氣渦輪火箭發(fā)動(dòng)機(jī)摻混燃燒研究進(jìn)展［J］.火箭推進(jìn)，2011，37(6):14-19.

［6］陳玉春，梁振欣，陳湘，等.固體推進(jìn)劑吸氣式渦輪發(fā)動(dòng)機(jī)控制規(guī)律與特性［J］.推進(jìn)技術(shù)，2009，30(3):332-336.

［7］梁振欣，陳玉春，黃興魯，等.SPATR發(fā)動(dòng)機(jī)循環(huán)參數(shù)選擇與特性分析［J］.固體火箭技術(shù)，2010，33(5):528-532.

［8］李江，王偉，劉洋，等.固體燃?xì)鉁u輪火箭發(fā)動(dòng)機(jī)摻混燃燒實(shí)驗(yàn)方法［J］.固體火箭技術(shù)，2013，36(2):170-174.

［9］Bossard J A，Christensen K L，F(xiàn)edun M H.Return of the solid fuel gas generator ATR［C］//23rd AIAA/SAE/ASME Joint Propulsion Conference.San Diego，California，1987.

［10］Calvo W C，Christensen K L，F(xiàn)edun M H.Solid fuel gas generator ATR［C］//22nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference.Huntsville，Alabama，1986.

［11］安德森約翰 D.計(jì)算流體力學(xué)基礎(chǔ)及其應(yīng)用［M］.吳頌平，譯.北京:機(jī)械工業(yè)出版社，2007.

［12］閻超.計(jì)算流體力學(xué)方法及應(yīng)用［M］.北京:北京航空航天大學(xué)出版社，2006.

［13］李平，李文龍，何國(guó)強(qiáng).空氣渦輪火箭發(fā)動(dòng)機(jī)內(nèi)外涵氣流摻混研究［J］.固體火箭技術(shù)，2012，35(1):47-52.

［14］Ross H.Design and testing of combustor for turbo-ramjet for UAV missile applications［D］.Montery，California:Naval Postgraduate School，2003.

［15］高嶺松，何國(guó)強(qiáng)，劉佩進(jìn)，等.固體火箭沖壓發(fā)動(dòng)機(jī)補(bǔ)燃室摻混與燃燒流場(chǎng)數(shù)值模擬［J］.固體火箭技術(shù)，2005，28(3):165-168.

［16］Haupt R L，Haupt S E.Practical genetic algorithms［M］.JoHn Wiley ＆ Son，Inc.，1998.

［17］田景文，高美娟.人工神經(jīng)網(wǎng)絡(luò)算法研究及應(yīng)用［M］.北京:北京理工大學(xué)出版社，2006.

［18］陳魁.試驗(yàn)設(shè)計(jì)與分析［M］.北京:清華大學(xué)出版社，2005.

［19］Engineous software inc.iSIGHT version 9.0 user＇s guide［Z］.Engineous Software，2004.