基于自適應(yīng)Kriging代理模型的葉型氣動(dòng)優(yōu)化設(shè)計(jì)

余 銳， 周代偉， 竺曉程， 杜朝輝

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

2.上海電氣電站設(shè)備有限公司上海汽輪機(jī)廠，上海200240）

葉型的氣動(dòng)設(shè)計(jì)是現(xiàn)代葉輪機(jī)械設(shè)計(jì)的核心技術(shù).合理高效的葉型是提高葉輪機(jī)械性能、降低葉輪機(jī)械運(yùn)行成本的前提條件之一.葉型的優(yōu)化設(shè)計(jì)往往要求在滿足各種性能結(jié)構(gòu)約束的前提下改善葉型的氣動(dòng)參數(shù).目前葉型的氣動(dòng)外形設(shè)計(jì)有反設(shè)計(jì)和優(yōu)化設(shè)計(jì)2種方法.反設(shè)計(jì)需要給出目標(biāo)壓力分布或表面速度分布，將目標(biāo)壓力分布轉(zhuǎn)化為葉型幾何外形是一個(gè)比較困難且對(duì)經(jīng)驗(yàn)有較高要求的問題.隨著優(yōu)化算法的發(fā)展和計(jì)算機(jī)水平的提高，基于隨機(jī)性尋優(yōu)的優(yōu)化設(shè)計(jì)方法被廣泛使用.雖然隨機(jī)性優(yōu)化算法具有良好的全局性、魯棒性及高并行性等特點(diǎn)，但多維大尋優(yōu)空間導(dǎo)致收斂速度慢是制約該算法應(yīng)用的主要障礙.建立高效的目標(biāo)評(píng)價(jià)方法、合理設(shè)計(jì)優(yōu)化參數(shù)空間和提高優(yōu)化算法的尋優(yōu)效率是應(yīng)用隨機(jī)性優(yōu)化算法的重要前提.雖然目前大多數(shù)透平型線已具有較高的熱力效率，但是作為組成葉片基本單元的型線，對(duì)其進(jìn)行優(yōu)化設(shè)計(jì)仍是目前進(jìn)行三維葉片設(shè)計(jì)的必要步驟.

基于梯度的優(yōu)化算法最早被運(yùn)用于葉片的氣動(dòng)數(shù)值優(yōu)化設(shè)計(jì)中.Burguburu［1］采用基于梯度的優(yōu)化算法耦合N－S求解器對(duì)3種葉輪機(jī)械葉片進(jìn)行了優(yōu)化設(shè)計(jì)，提高了葉片的效率.然而由于基于梯度的優(yōu)化算法易陷入局部最優(yōu)，且難以適應(yīng)葉片氣動(dòng)優(yōu)化設(shè)計(jì)多模態(tài)的特點(diǎn)，因而一些具有較好全局尋優(yōu)能力的啟發(fā)式算法逐漸廣泛應(yīng)用于優(yōu)化設(shè)計(jì)中［2］.Oyama［3］采用遺傳算法對(duì) NASA Rotor67壓氣機(jī)葉片進(jìn)行優(yōu)化設(shè)計(jì)，取得了較好的效果.?KSüZ［4］介紹了一種用于三維燃?xì)廨啓C(jī)葉片的氣動(dòng)優(yōu)化設(shè)計(jì)，采用遺傳算法分別對(duì)葉片的氣動(dòng)效率和扭矩進(jìn)行單目標(biāo)優(yōu)化，結(jié)果表明，優(yōu)化后葉片的氣動(dòng)效率從83.9%提高到85.9%，扭矩增大了7.6%，完成了對(duì)原型葉片的改進(jìn).全局尋優(yōu)的算法需要大量CFD計(jì)算，時(shí)間耗費(fèi)巨大，不能滿足實(shí)際的氣動(dòng)設(shè)計(jì)和工程應(yīng)用.在設(shè)計(jì)變量與優(yōu)化目標(biāo)間構(gòu)造近似模型（又稱代理模型）來替代CFD計(jì)算進(jìn)行尋優(yōu)可以很好地解決上述問題，因此，該方法得到了越來越多的應(yīng)用［5］.薛亮等［6］采用二次多項(xiàng)式近似模型及改進(jìn)遺傳算法，省時(shí)高效地完成了對(duì)NASA Rotor67動(dòng)葉片葉型的優(yōu)化設(shè)計(jì)，使得葉片絕熱效率提高了1.1%.高行山等［7］在渦輪葉片氣動(dòng)優(yōu)化時(shí)引入近似技術(shù)，加速了整個(gè)尋優(yōu)過程，通過對(duì)某渦輪葉片的優(yōu)化分析驗(yàn)證了該方法的可行性.Mengistu等［8］采用遺傳算法構(gòu)造優(yōu)化目標(biāo)與優(yōu)化變量間的近似模型，以最大化絕熱效率、壓比和最小化總壓損失為優(yōu)化目標(biāo)，對(duì)亞音速渦輪葉柵二維流動(dòng)情況進(jìn)行優(yōu)化改進(jìn).最后通過4組方案分析了設(shè)計(jì)參數(shù)對(duì)葉片氣動(dòng)性能的影響，得到了更優(yōu)的葉型.王紅濤［9］對(duì)比分析了基于自適應(yīng)Kriging代理模型的優(yōu)化算法與小生境微種群遺傳算法，發(fā)現(xiàn)基于自適應(yīng)Kriging代理模型的引入減少了迭代所需調(diào)用CFD計(jì)算的次數(shù)，有效地提高了優(yōu)化過程的效率.

筆者采用基于自適應(yīng)Kriging代理模型的近似代理技術(shù)建立具有全局尋優(yōu)能力的優(yōu)化算法，通過對(duì)葉型進(jìn)行參數(shù)化描述，代入RANS方程求解得到葉型氣動(dòng)性能參數(shù)，建立葉型氣動(dòng)優(yōu)化設(shè)計(jì)方法，并以壓力恢復(fù)系數(shù)為目標(biāo)函數(shù)，使用該方法對(duì)某設(shè)計(jì)工況下的葉型進(jìn)行了優(yōu)化設(shè)計(jì)，驗(yàn)證了整個(gè)優(yōu)化系統(tǒng)的可靠性，提高了葉型的性能.

1 數(shù)值優(yōu)化方法

一般對(duì)于帶有約束條件的優(yōu)化問題可采用如下表達(dá)形式

式中：X為設(shè)計(jì)變量；fopt（X）為目標(biāo)函數(shù)的最優(yōu)值；X＝［x1，x2，…，xn］T，代表在n維空間Rn中的點(diǎn)，即優(yōu)化問題中的1組方案；Ω為設(shè)計(jì)變量的定義域；gi（X）和hj（X）為約束函數(shù).

式（2）中的優(yōu)化問題需要滿足2種約束條件：不等式約束和等式約束.按照規(guī)范化形式，一般都將優(yōu)化問題歸結(jié)為目標(biāo)函數(shù)的于極小值問題，對(duì)于目標(biāo)函數(shù)極大化問題可以通過轉(zhuǎn)換為對(duì)應(yīng)的極小化問題求解.

圖1給出了基于自適應(yīng)Kriging代理模型的葉型氣動(dòng)優(yōu)化設(shè)計(jì)流程示意圖.由圖1可知，該系統(tǒng)以調(diào)節(jié)葉型形狀的控制參數(shù)作為設(shè)計(jì)變量，通過試驗(yàn)設(shè)計(jì)（DOE）獲得初始的設(shè)計(jì)變量樣本集.并在完成對(duì)葉型網(wǎng)格的自動(dòng)生成和流場(chǎng)的求解分析后，獲取目標(biāo)函數(shù)信息（即響應(yīng)樣本集）.進(jìn)而采用自適應(yīng)Kriging代理模型構(gòu)造變量和響應(yīng)間的近似模型，完成校正點(diǎn)的尋找和評(píng)估.若不收斂，將校正點(diǎn)和對(duì)應(yīng)響應(yīng)值加入樣本庫，不斷更新代理模型直至收斂，獲得目標(biāo)函數(shù)最優(yōu)的葉型.該系統(tǒng)將試驗(yàn)設(shè)計(jì)、CFD數(shù)值模擬和優(yōu)化設(shè)計(jì)方法相結(jié)合，自動(dòng)完成各模塊間的數(shù)據(jù)傳遞，減少了人工干預(yù)且極大地提高了葉型開發(fā)效率.

圖1 基于自適應(yīng)Kriging代理模型的葉型氣動(dòng)優(yōu)化設(shè)計(jì)流程Fig.1 Flowchart of airfoil aerodynamic optimization design based on adaptive Kriging surrogate model

1.1 優(yōu)化算法

全局優(yōu)化算法（如遺傳算法等）雖然可以使優(yōu)化問題避免陷入局部最優(yōu)，但由于耗時(shí)巨大，難以直接應(yīng)用于實(shí)際工程優(yōu)化設(shè)計(jì)中.因而采用代理模型代替數(shù)值模擬手段是一種行之有效的方法，可以大大減少優(yōu)化過程的計(jì)算量、提高優(yōu)化效率及優(yōu)化算法的性能和組成更好的優(yōu)化策略.

常用的代理模型主要有多項(xiàng)式響應(yīng)面模型、徑向基函數(shù)模型和Kriging模型等.筆者對(duì)葉型優(yōu)化時(shí)選取的算法是一種基于自適應(yīng)Kriging代理模型的自適應(yīng)序貫優(yōu)化算法［10］.由于Kriging模型在訓(xùn)練樣本點(diǎn)處滿足無偏估計(jì)，具有良好的高度非線性近似能力，因此最適合在葉型優(yōu)化等高度非線性優(yōu)化問題中作為代理模型使用.

Kriging模型假設(shè)樣本點(diǎn)x（維數(shù)為m）與對(duì)應(yīng)響應(yīng)值y之間的關(guān)系為

式中：f（）x為確定性部分，可以用一個(gè)常數(shù)表示，是對(duì)設(shè)計(jì)空間的全局近似；z（x）為一隨機(jī)過程，表示對(duì)全局近似的背離.

該隨機(jī)過程與全局模型的局部偏差是由相關(guān)函數(shù)反映的.z（x）具有期望為0，方差為σ2的統(tǒng)計(jì)特征.方差σ2反映了試驗(yàn)的精度，一般σ2越大，表示試驗(yàn)誤差越大.隨機(jī)過程z（x）的協(xié)方差為：

式中：R為相關(guān)函數(shù)矩陣；ns為樣本點(diǎn)個(gè)數(shù)；R （xi，xj）為點(diǎn)（xi，xj）間的相關(guān)函數(shù)，此相關(guān)函數(shù)為高斯函數(shù).

相關(guān)函數(shù)的形式如下：

式中：θk為各向異性參數(shù).

在預(yù)測(cè)點(diǎn)x0處的預(yù)測(cè)響應(yīng)值和預(yù)測(cè)標(biāo)準(zhǔn)差分別為

相關(guān)函數(shù)R （xi，xj）及和均依賴于各向異性參數(shù)θ，因而構(gòu)建Kriging模型的關(guān)鍵在于θ的求解，通常由極大似然估計(jì)給出，即在θ＞0的情況下使式（8）取得極大值.通過求解得到各向異性參數(shù)θ后，即完成了整個(gè)Kriging代理模型的構(gòu)建.

1.2 參數(shù)化葉型

對(duì)葉型進(jìn)行參數(shù)化表達(dá)是實(shí)現(xiàn)自動(dòng)優(yōu)化設(shè)計(jì)的前提，一般二維葉型的參數(shù)化主要有2種方法：中弧線加厚度分布和采用構(gòu)造曲線直接對(duì)葉型壓力面、吸力面型線進(jìn)行參數(shù)化.目前常用的構(gòu)造曲線主要有3種：Bézier曲線、均勻有理B樣條曲線和非均勻有理B樣條即NURBS曲線［11］.其中Bézier曲線由于形式簡(jiǎn)單、易操作，在參數(shù)化中得到廣泛應(yīng)用.

采用2段Bézier曲線分別對(duì)某渦輪葉型壓力面和吸力面型線進(jìn)行參數(shù)化，共有16個(gè)控制點(diǎn)控制葉型形狀.其中位于原點(diǎn)處的控制點(diǎn)為2段Bézier曲線的公共點(diǎn)，吸力面型線由自公共點(diǎn)順時(shí)針連接的10個(gè)控制點(diǎn)獲得，而壓力面型線由自公共點(diǎn)逆時(shí)針連接的7個(gè)控制點(diǎn)獲得.圖2為葉型型線及對(duì)應(yīng)控制點(diǎn)連線示意圖，其中x和y分別表示葉型的軸向和周向位置.對(duì)葉型進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)時(shí)，通過改變控制點(diǎn)位置即可生成新的葉型，直到獲得設(shè)計(jì)變量空間下目標(biāo)函數(shù)最優(yōu)時(shí)的葉型，完成葉型優(yōu)化改進(jìn).

圖2 葉型型線及對(duì)應(yīng)控制點(diǎn)連線Fig.2 Profile of the airfoil and polyline of airfoil control points

1.3 數(shù)值方法

在網(wǎng)格生成時(shí)選用HOH型拓?fù)浣Y(jié)構(gòu)：進(jìn)口段為H型結(jié)構(gòu)化網(wǎng)格，葉型處為O型結(jié)構(gòu)化網(wǎng)格，考慮到實(shí)際流動(dòng)情況，出口段選擇I型網(wǎng)格.對(duì)近壁面處網(wǎng)格進(jìn)行加密，使葉片壁面第1層網(wǎng)格滿足y＋＜10.圖3給出了葉型前緣和尾緣處的網(wǎng)格.

圖3 葉型前緣和尾緣處的網(wǎng)格Fig.3 Mesh for the leading and trailing edge of airfoil

對(duì)要進(jìn)行氣動(dòng)優(yōu)化研究的渦輪葉型構(gòu)造平面葉柵，葉柵柵距為50mm.數(shù)值計(jì)算時(shí)給定進(jìn)出口邊界條件，其中進(jìn)口處給定平均總壓為2×106Pa，進(jìn)口總溫為573K，進(jìn)口來流氣流角為5°；在出口處給定平均靜壓，使得出口靜壓與進(jìn)口總壓比為0.896.湍流模型選擇一方程的Spalart－Allmaras湍流模型（即S－A湍流模型），空間離散采用二階精度的迎風(fēng)格式，時(shí)間項(xiàng)推進(jìn)采用四階Runge－Kutta法加速迭代求解.最大迭代步數(shù)設(shè)為5 000，收斂標(biāo)準(zhǔn)設(shè)為全局殘差降到10－6以下.

為驗(yàn)證本文數(shù)值方法的可靠性，選取文獻(xiàn)［12］中出口氣流角為67°的直葉片進(jìn)行校核.該文獻(xiàn)采用試驗(yàn)手段和理論計(jì)算對(duì)不同出口氣流角的葉柵進(jìn)行氣動(dòng)性能的研究，給出了詳細(xì)的葉片幾何數(shù)據(jù)和試驗(yàn)數(shù)據(jù).根據(jù)上述數(shù)據(jù)生成葉柵網(wǎng)格，采用S－A湍流模型和二階精度的迎風(fēng)格式在設(shè)計(jì)工況下進(jìn)行數(shù)值模擬.圖4給出了中徑處葉片表面臨界速度比沿相對(duì)軸向弦長(zhǎng)的分布對(duì)比.由圖4可以看出，數(shù)值計(jì)算所得結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好，從而驗(yàn)證了數(shù)值計(jì)算方法的可靠性.

圖4 中徑處葉片表面臨界速度比分布對(duì)比Fig.4 Comparison of surface critical velocity ratio at mid－span height between experimental and numerical data

2 渦輪葉型優(yōu)化設(shè)計(jì)

采用2組Bézier控制點(diǎn)控制型線變化，其中O點(diǎn)為壓力面型線和吸力面型線的公共控制點(diǎn)，如圖5所示.在生產(chǎn)實(shí)踐中，常通過微調(diào)葉片的頭部型線來滿足不同設(shè)計(jì)要求下進(jìn)出口流動(dòng)角的變化，為此需要研究工程設(shè)計(jì)中對(duì)頭部型線的微調(diào)是否合理及是否存在最優(yōu)頭部型線.以對(duì)葉型頭部控制點(diǎn)的優(yōu)化為例，采用上述氣動(dòng)優(yōu)化設(shè)計(jì)方法，研究其對(duì)葉型氣動(dòng)性能的優(yōu)化設(shè)計(jì).設(shè)定原始葉型控制點(diǎn)A和B到O點(diǎn)的距離分別為R1和R2，控制點(diǎn)變化后對(duì)應(yīng)的距離為r1和r2.優(yōu)化時(shí)為保證葉型幾何進(jìn)氣角不變，只改變OA和OB2段控制線的長(zhǎng)度而不改變其角度，即以系數(shù)矩陣L作為設(shè)計(jì)變量，以壓力恢復(fù)系數(shù)ξ最大為優(yōu)化目標(biāo)進(jìn)行葉型的優(yōu)化設(shè)計(jì).

圖5 葉型前緣控制點(diǎn)Fig.5 Control points at leading edge of profile

系數(shù)矩陣L與壓力恢復(fù)系數(shù)ξ的定義如下式中：p01和p02分別為平面葉柵進(jìn)口和出口的總壓；p2為出口靜壓.

渦輪葉型氣動(dòng)優(yōu)化問題可以用式（11）來表示

利用試驗(yàn)設(shè)計(jì)方法和U＊10（108）均勻設(shè)計(jì)表生成10個(gè)初始樣本點(diǎn)，使樣本點(diǎn)在試驗(yàn)空間里均勻分布，并構(gòu)造初始化代理模型，進(jìn)而對(duì)優(yōu)化目標(biāo)進(jìn)行尋優(yōu).圖6給出了整個(gè)優(yōu)化過程中響應(yīng)值隨迭代步數(shù)的變化.

圖6 響應(yīng)值隨迭代步數(shù)的變化Fig.6 Response of pressure recovery coefficient to iteration step

設(shè)計(jì)工況下原始葉型的壓力恢復(fù)系數(shù)ξ0＝0.973 76，通過優(yōu)化迭代后得到當(dāng)設(shè)計(jì)變量L＝［1.800 000，1.032 544］時(shí)，總壓恢復(fù)系數(shù)達(dá)到最大，即ξ1＝0.974 40，壓力恢復(fù)系數(shù)有所增大.優(yōu)化結(jié)果表明，存在最優(yōu)頭部型線，但壓力恢復(fù)系數(shù)增大的程度有限；同時(shí)也說明了實(shí)際工程設(shè)計(jì)中所采用的頭部型線微調(diào)方法的合理性.圖7給出了優(yōu)化前后控制點(diǎn)及型線的對(duì)比.從圖7可以看出，葉型前緣處控制線段長(zhǎng)度優(yōu)化的主要變化集中在吸力面型線上，壓力面型線整體上保持不變.

為進(jìn)一步分析葉型前緣型線的變化對(duì)流場(chǎng)的影響，圖8給出了原始葉型和優(yōu)化葉型表面壓力系數(shù)Cp沿軸向的分布，其中表面壓力系數(shù)Cp定義如下

式中：p為當(dāng)?shù)仂o壓；p0為葉柵進(jìn)口平均靜壓；p＊0為葉柵進(jìn)口平均總壓.

圖7 優(yōu)化前后控制點(diǎn)及型線的對(duì)比Fig.7 Comparison between original and optimized profile

圖8 原始葉型和優(yōu)化葉型表面壓力系數(shù)沿軸向的分布Fig.8 Comparison of pressure coefficient between original and optimized airfoil

從圖8可以看出，優(yōu)化前后葉型壓力面上的表面壓力系數(shù)分布規(guī)律基本一致；而在吸力面上存在明顯的差異，優(yōu)化葉型靠近前緣吸力面處的表面壓力系數(shù)較小，氣流沿吸力面流動(dòng)時(shí)壓力變化較平緩.葉型前緣變化引起的流動(dòng)變化使得葉型型面損失減小，壓力恢復(fù)系數(shù)有所增大，達(dá)到了優(yōu)化葉型、提高葉柵性能的目的.

3 結(jié) 論

基于自適應(yīng)Kriging代理模型，設(shè)計(jì)開發(fā)了一套集參數(shù)化建模、優(yōu)化算法和數(shù)值求解于一體的自動(dòng)葉型氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng).利用該系統(tǒng)對(duì)某渦輪葉型頭部型線進(jìn)行了優(yōu)化設(shè)計(jì)，使原始葉型的壓力恢復(fù)系數(shù)進(jìn)一步增大，結(jié)果表明存在最優(yōu)頭部型線，說明了工程實(shí)踐中對(duì)頭部型線調(diào)整的合理性.經(jīng)過優(yōu)化改善了葉型的性能，同時(shí)也驗(yàn)證了整個(gè)優(yōu)化系統(tǒng)的可靠性，為葉片的全三維優(yōu)化設(shè)計(jì)、提高葉片及級(jí)性能的研究奠定了基礎(chǔ).

［1］BURGUBURU S，TOUSSAINT C，BONHOMME C.Numerical optimization of turbomachinery bladings［J］.Journal of Turbomachinery，2004，126（1）：91－100.

［2］ZHOU Fangzhen，F(xiàn)ENG Guotai，JIANG Hongde.The development of highly loaded turbine rotating blades by using 3doptimization design method of turbomachinery blades based on artificial neural network＆genetic algorithm［J］.Chinese Journal of Aeronautics，2003，16（4）：198－202.

［3］OYAMA A，LIOU M S，OBAYASHI S.Transonic axial flow blade optimization：evolutionary algorithms／three dimensional Navier－Stokes solver［J］.Journal of Propulsion and Power，2004，20（4）：612－619.

［4］?KSüZ?，AKMANDOR˙I.Turbine blade shape aerodynamic design using artificial intelligence［C］／／ASME Turbo Expo.Reno－Tahoe，Nevada：ASME，2005.

［5］王紅濤，竺曉程，杜朝輝.自適應(yīng)Kriging近似模型及其在二維擴(kuò)壓器優(yōu)化設(shè)計(jì)中的應(yīng)用［J］.計(jì)算力學(xué)學(xué)報(bào)，2011，28（1）：15－19，24.WANG Hongtao，ZHU Xiaocheng，DU Zhaohui.Applications of adaptive Kriging approximation model in two dimensional diffuser aerodynamic optimization design［J］.Chinese Journal of Computational Mechanics，2011，28（1）：15－19，24.

［6］薛亮，韓萬金.基于遺傳算法與近似模型的全局氣動(dòng)優(yōu)化方法［J］.推進(jìn)技術(shù)，2008，29（3）：21－26.XUE Liang，HAN Wanjin.Global aerodynamic optimization method using genetic algorithms and surrogate model［J］.Journal of Propulsion Technology，2008，29（3）：21－26.

［7］高行山，韓永志，張娟，等.基于近似技術(shù)的渦輪葉片氣動(dòng)優(yōu)化設(shè)計(jì)［J］.計(jì)算力學(xué)學(xué)報(bào)，2008，25（6）：874－877.GAO Hangshan，HAN Yongzhi，ZHANG Juan，et al.Blade aerodynamic optimization design based on approximation method［J］.Chinese Journal of Computational Mechanics，2008，25（6）：874－877.

［8］MENGISTU T，GHALY W，MANSOUR T．Aerodynamic shape optimization of turbine blades using a design－parameter－based shape representation［C］／／ASME Turbo Expo．Montreal，Canada：ASME，2007．

［9］王紅濤．汽輪機(jī)低壓排汽系統(tǒng)內(nèi)部流動(dòng)及其氣動(dòng)優(yōu)化設(shè)計(jì)研究［D］．上海：上海交通大學(xué)，2011．

［10］王紅濤，竺曉程，杜朝輝．基于Kriging代理模型的自適應(yīng)序貫優(yōu)化方法［J］．計(jì)算機(jī)工程與應(yīng)用，2009，45（26）：193－195．WANG Hongtao，ZHU Xiaocheng，DU Zhaohui．A－daptive sequential optimization algorithm based on Kriging surrogate model［J］．Computer Engineering and Applications，2009，45（26）：193－195．

［11］樂英，王璋奇，韓慶瑤．基于非均勻有理B樣條曲面重構(gòu)汽輪機(jī)扭曲葉片［J］．動(dòng)力工程學(xué)報(bào)，2010，30（12）：926－931．YUE Ying，WANG Zhangqi，HAN Qingyao．Reconstruction of turbine blade twisted based on NURBS surface［J］．Journal of Chinese Society of Power Engineering，2010，30（12）：926－931．

［12］SCHWAB J R．Aerodynamic performance of high turning core turbine vanes in a two－dimensional cascade［C］／／AIAA／SAE／ASME 18th Joint Propulsion Conference．Cleveland，Ohio：AIAA，1982．