基于CFD的多級推力固體火箭發(fā)動機軸對稱噴管型面優(yōu)化與高精度性能預(yù)估①

陳 偉，梁國柱

(北京航空航天大學宇航學院，北京 100191)

0 引言

伴隨著導(dǎo)彈武器系統(tǒng)對固體火箭發(fā)動機性能越來越高的要求，噴管的工作效率對固體火箭發(fā)動機性能的影響也就占有更加重要的位置。因此，如何高水平地優(yōu)化噴管型面，特別是在現(xiàn)有噴管的基礎(chǔ)上，如何進一步深入挖掘噴管的工作潛能，就成為固體火箭發(fā)動機設(shè)計的重要課題。

噴管型面的優(yōu)化設(shè)計問題，首先要解決的問題是如何高效準確地預(yù)估噴管型面的性能和效率。早期的優(yōu)化設(shè)計受計算能力的限制，大多采用計算量小、且模型相對簡單(無粘流)的特征線方法［1］，而隨著計算流體力學技術(shù)和并行計算技術(shù)的發(fā)展，現(xiàn)在越來越多的噴管型面優(yōu)化研究是基于更為精確的直接數(shù)值求解N-S方程的 CFD仿真方法［2－3］。但由于 CFD方法復(fù)雜、計算量大計算穩(wěn)定性差等問題，若不采取有效措施和策略，很難將CFD仿真方法和噴管的優(yōu)化設(shè)計直接結(jié)合在一起。因此，如何能夠高效穩(wěn)健地基于CFD仿真方法來對噴管內(nèi)型面進行優(yōu)化設(shè)計，是噴管的優(yōu)化設(shè)計中較為重要的課題。

本文將以較為常見的多級推力的固體火箭發(fā)動機為例，以CFD仿真為基礎(chǔ)，探索合適的優(yōu)化方法和策略，達到高效穩(wěn)健和準確地優(yōu)化設(shè)計噴管型面的目的，同時實現(xiàn)噴管的高精度性能預(yù)估。

1 優(yōu)化對象和策略

1．1 優(yōu)化問題分析

固體火箭發(fā)動機噴管內(nèi)的流動由于摩擦、兩相流等原因，使得噴管流動過程中存在各種損失，噴管實際比沖小于理論比沖。通過優(yōu)化噴管型面，能降低噴管流動損失，提高發(fā)動機性能［4］。

多級推力固體火箭發(fā)動機相對于單級推力固體火箭發(fā)動機而言，其內(nèi)彈道曲線和推力曲線呈階梯狀。如圖1所示的雙推力發(fā)動機的內(nèi)彈道曲線，其中每個相對平穩(wěn)的工作時間段稱為工作段(圖1中I和II)。多級推力的固體火箭發(fā)動機在導(dǎo)彈武器系統(tǒng)中是一種較為常見的發(fā)動機類型。例如，空空導(dǎo)彈的研制中，就廣泛采用單室雙推力的固體火箭發(fā)動機，因為雙推力對應(yīng)的導(dǎo)彈平均速度與末端速度都要比單推力大，性能更好［5］。由于燃燒室壓力在不同的工作段各不相同，因此噴管在各個階段的效率也各不相同，即使在某一個工作段噴管達到最優(yōu)效率，往往其他工作段的效率就相對很低，整個工作段的效率更不可能達到最優(yōu)。因此，需要兼顧各個階段的噴管效率，才能使得噴管工作過程達到最優(yōu)性能。而對噴管工作性能的評價，可采用CFD仿真方法來計算噴管內(nèi)流場，從而通過分析流場的物理參數(shù)，來預(yù)測噴管的實際工作性能。

圖1 典型的雙推力固體火箭發(fā)動機的內(nèi)彈道曲線Fig．1 Typical interior ballistics curve of dual-stage-thrust SRM

解決多級推力固體火箭發(fā)動機軸對稱噴管的優(yōu)化問題，首先要提取噴管型面的設(shè)計變量，建立CFD仿真模型;然后，通過CFD仿真方法預(yù)測多級推力固體火箭發(fā)動機各個不同工作階段噴管的工作性能，綜合各個工作段的噴管性能，預(yù)測噴管的總體工作性能;最后，通過優(yōu)化算法優(yōu)化型面的設(shè)計變量，使得噴管的性能達到最優(yōu)。

1．2 仿真和優(yōu)化策略

無論是CFD仿真計算，還是數(shù)值優(yōu)化計算，都是通過迭代計算的方法來得到計算結(jié)果的，而基于CFD的優(yōu)化就相當于嵌套迭代，其計算量相當于單獨計算量的二階放大。如果不采取措施，其計算量會達到難以承受的地步。另一方面，CFD仿真計算和數(shù)值優(yōu)化計算都存在收斂穩(wěn)健性問題，在很多情況下會出現(xiàn)迭代無法收斂的情況，而基于CFD仿真的優(yōu)化，其收斂穩(wěn)健性會變得更為脆弱。因此，必須采用合理的措施和策略，在保證計算精度和結(jié)果可靠性的前提下，減少計算規(guī)模，提高計算效率，且要保證計算的穩(wěn)健性，提高計算的可靠性。根據(jù)本研究優(yōu)化對象的特點，本文從CFD仿真模型、仿真方法、優(yōu)化流程3個方面，來研究提高計算效率和計算穩(wěn)健性的優(yōu)化和仿真策略。

首先，從以下2個方面來簡化CFD仿真模型，減少計算規(guī)模。

(1)將內(nèi)彈道曲線分解，燃燒室壓力在各個工作段上的變化較為平緩，可近似認為在各個工作段內(nèi)為定常流。因此，采用定常流模型來進行CFD仿真，每個工作段的平均壓強作為典型噴管入口壓強進行仿真計算。

(2)固體火箭發(fā)動機噴管內(nèi)流場多為氣固兩相流，但經(jīng)過實踐表明，采用兩相流模型進行仿真的計算時間，往往是在相同條件下采用兩相平衡流模型的計算時間的2～3倍，且兩相流模型的計算穩(wěn)健性要差于兩相平衡流模型。因此，采用結(jié)合兩相平衡流計算模型和兩相流計算模型的方式來保證仿真精度，同時減少計算量和提高計算穩(wěn)健性:優(yōu)化前期，采用兩相平衡流模型對噴管型面進行優(yōu)化;然后，在兩相平衡流的優(yōu)化型面基礎(chǔ)上，利用兩相流模型對其進行修正優(yōu)化。

(3)傳統(tǒng)CFD仿真的一般流程是根據(jù)仿真對象建立2維或3維幾何模型，根據(jù)問題劃分網(wǎng)格，選擇和配置CFD仿真模型，然后進行仿真迭代計算。這種方法并不適合于優(yōu)化中的仿真計算，因為每次重新劃分網(wǎng)格很難做到仿真的自動化，也無法保證網(wǎng)格的質(zhì)量。考慮到本文研究的軸對稱噴管是在其基本結(jié)構(gòu)和工作狀態(tài)大體已經(jīng)確定的情況下的進一步挖潛優(yōu)化，在各種幾何約束的嚴格限制下，噴管型面尺寸參數(shù)可變的幅度已經(jīng)很小。因此，本文采用了動態(tài)網(wǎng)格技術(shù)，在基本不改變網(wǎng)格規(guī)模的前提下，通過調(diào)整網(wǎng)格分布來適應(yīng)新型面的方法，來處理每次的網(wǎng)格生成，使得每次新型面的仿真可采用前一次型面的流場數(shù)據(jù)作為仿真的初始值。由于每次優(yōu)化中產(chǎn)生型面變化不大，流場改變很小。這樣極大地減少了每次仿真的計算量，也大大提高了仿真計算穩(wěn)健性。具體方法見第5章。

(4)由于各個工作段的仿真計算是相互獨立的，因此可將各個工作段的仿真進行并行計算，提高計算效率。另外，在優(yōu)化算法的選取上，各個優(yōu)化算法的特點各不相同，有的收斂速度快，迭代穩(wěn)健性差，而有的收斂速度慢，但迭代穩(wěn)健性較好。因此，在本研究中，采用多個優(yōu)化算法相結(jié)合的辦法，在保證優(yōu)化的穩(wěn)健性基礎(chǔ)上，盡量提高優(yōu)化的收斂速度，減少計算量。

下面分別對仿真和優(yōu)化模型進行具體闡述。

2 CFD仿真模型

2．1 基本假設(shè)

依據(jù)前文對CFD仿真模型的分析，做出如下假設(shè):

(1)噴管內(nèi)流場為二維軸對稱定常流;

(2)燃氣中的氣相為完全氣體，服從完全氣體狀態(tài)方程;

(3)忽略熱輻射作用和化學反應(yīng)，近似認為噴管內(nèi)流動過程絕熱;

(4)顆粒為均勻球體，顆粒是離散的，不考慮顆粒相的燃燒、蒸發(fā)、破碎及化學反應(yīng)。不考慮由于顆粒相加速度產(chǎn)生的虛擬質(zhì)量力，不考慮外部體積力和升力作用，不考慮顆粒相的體積分數(shù)變化湍流脈動產(chǎn)生的壓力。

2．2 控制方程

采用Euler形式的二維軸對稱有粘兩相定常流的控制方程［6］:

連續(xù)方程為

pqq相的體積分數(shù);ρp和ρq分別為第p相和第q相的密度。不考慮兩相之間的傳質(zhì)現(xiàn)象。

動量守恒方程為

其中，Θs為顆粒溫度;ess為顆粒碰撞恢復(fù)系數(shù);g0，ss為顆粒徑向分布函數(shù)。而對于氣相，pq=0。

其中，s表示顆粒相參數(shù);l表示氣相參數(shù);τs表示顆粒相微粒弛豫時間，定義為

阻力系數(shù)CD為

能量方程為

式中hq為第q相的比焓為第q相熱流矢量。

注意到，若將式(1)～式(7)中固相的相應(yīng)參數(shù)設(shè)為0，即可得到兩相平衡流CFD仿真模型的控制方程。

采用較成熟的標準k-ε兩方程湍流模型和重點區(qū)域加密的結(jié)構(gòu)網(wǎng)格，以提高計算的準確性和穩(wěn)健性。

2．3 仿真算法和邊界條件

兩相平衡流采用基于密度的隱式耦合算法，并選擇二階迎風的Roe-FDS格式作為數(shù)值離散格式。兩相流采用多相耦合算法，即同時求解所有的方程來得到相速度和相間共享壓力的修正值［9］，離散格式均采用一階迎風格式。

入口邊界條件為壓力入口，其壓力為各個階段壓力的平均值。對于兩相平衡流，入口總溫為燃燒室燃氣總溫。對于氣固兩相流，氣相壓力和顆粒相均采用兩相平衡流的燃燒室燃氣總溫作為入口總溫，而顆粒相的入口體積分數(shù)由推進劑顆粒相質(zhì)量分數(shù)和氣相與顆粒相的入口密度換算得到。

出口邊界條件為壓力出口，靜壓為其發(fā)動機工作高度的大氣壓，靜溫為室溫。但由于噴管出口處為超聲速氣流，噴管出口的壓力和溫度由內(nèi)部氣流推斷，而非邊界條件上的壓力和溫度。當出口的局部壓力低于邊界條件的靜壓時，會發(fā)生回流現(xiàn)象，這時回流的壓強和溫度為邊界條件的靜壓和靜溫。

另外，兩相流顆粒大小采用質(zhì)量加權(quán)平均經(jīng)驗公式和體積-表面積平均半徑經(jīng)驗公式的平均值［10］:

3 噴管性能預(yù)估

在由CFD仿真計算得到的單個工作段的噴管內(nèi)流場的基礎(chǔ)上，通過對噴管的推力和質(zhì)量流量進行積分，可計算出噴管工作的實際比沖，即

式中 Re為噴管出口半徑;F為推力;為質(zhì)量流量;Isi為第i個推力區(qū)間的噴管實際比沖。

因推進劑質(zhì)量mi是給定的，則根據(jù)前面的假設(shè)，可得到發(fā)動機的實際總沖為

引入熱力計算得到的噴管入口的燃氣參數(shù)，可直接計算得到各個工作段不計噴管損失的理想比沖，并進一步計算得到理想總沖，即

由于噴管入口的燃氣參數(shù)是給定的，因此計算噴管的總工作效率為實際比沖除以理想比沖，即

通過以上步驟計算得到噴管總工作效率η和實際總沖I指標，就可對多級推力的固體火箭發(fā)動機的噴管性能做出精確預(yù)估和評價。

4 優(yōu)化模型

4．1 優(yōu)化變量

噴管型面作為優(yōu)化的對象，需采用一種描述噴管型面的方法來提取優(yōu)化變量，本文采用分段線的描述法，即提取型面曲線上的特征變量。以圖2為例，這是一個典型拋物線噴管，將其分為5段曲線，從中可提取7個相互獨立的控制變量(其中，噴管的喉部直徑由發(fā)動機確定的，而噴管入口由燃燒室結(jié)構(gòu)設(shè)計所確定)，型面的其余變量包括噴管出口半徑，可由這7個變量推導(dǎo)得出，如表1所示。

4．2 數(shù)學模型

優(yōu)化的數(shù)學模型主要包含3個部分，即優(yōu)化變量、約束條件和目標函數(shù)。

表1 噴管的控制變量Table 1 Control variable of nozzle

圖2 拋物線噴管控制變量Fig．2 Control variable of parabola nozzle

優(yōu)化變量:本文的優(yōu)化對象是噴管型面，優(yōu)化變量為4．1節(jié)所描述型面的控制參數(shù)x，即

約束條件:噴管型面的主要約束包括結(jié)構(gòu)約束和形狀約束，結(jié)構(gòu)約束包括發(fā)動機其他結(jié)構(gòu)部分以及導(dǎo)彈/火箭結(jié)構(gòu)的限制和約束，而形狀約束主要是指噴管形狀的幾何合理性的約束，兩者均為不等式約束，可表示為

目標函數(shù):選擇由式(10)得到性能預(yù)估指標實際總沖I作為目標函數(shù):

考慮到優(yōu)化計算中一般以最小化目標，確定優(yōu)化模型為

4．3 優(yōu)化算法

本研究針對的噴管優(yōu)化問題是非常典型的非線性有約束的線性規(guī)劃問題。按照前述的優(yōu)化算法的策略，本文選擇了內(nèi)點法(Internal Point)［11］和序列二次規(guī)劃(SQP)［12］2種特點互補的算法完成該優(yōu)化計算。內(nèi)點法穩(wěn)健性高，收斂性好，但收斂速度慢，而SQP算法收斂速度高，但易發(fā)散。優(yōu)化計算前期，采用內(nèi)點法，保證其收斂穩(wěn)定性，在接近最優(yōu)點時，切換成SQP算法，提高收斂速度。利用2種算法相結(jié)合的辦法，提高收斂速度，保證了優(yōu)化迭代的穩(wěn)健性。

5 計算流程實現(xiàn)

5．1 動態(tài)網(wǎng)格生成

根據(jù)前文所述的策略，通過采用動態(tài)網(wǎng)格的方法，避免了每次重新對新型面進行CFD仿真，從而減少計算量、提高計算穩(wěn)健度。

本文采用Thompson提出的橢圓型方程的結(jié)構(gòu)網(wǎng)格生成方法［13］，在迭代初期，確定加密位置和網(wǎng)格總數(shù)，利用橢圓型方程結(jié)構(gòu)網(wǎng)格生成法，根據(jù)新型面參數(shù)調(diào)整網(wǎng)格點分布，既保證了網(wǎng)格的相似性，也保證了網(wǎng)格適用于新型面。流程圖如圖3所示。

圖3 動態(tài)網(wǎng)格流程圖Fig．3 Flow chart of dynamic grid generation

每次優(yōu)化迭代產(chǎn)生新的設(shè)計變量，通過新的設(shè)計變量，生成型面曲線，若是第一次生成網(wǎng)格，則由Thompson網(wǎng)格生成法生成新的網(wǎng)格進行CFD仿真計算。若不是，則通過Thompson網(wǎng)格生成法，在不改變網(wǎng)格總數(shù)和分布要求的基礎(chǔ)上，按照新型面參數(shù)生成網(wǎng)格，將上一個網(wǎng)格的模型參數(shù)、算法設(shè)置、流場結(jié)果等映射到新生成的網(wǎng)格，使其仍能使用上一次仿真的流場計算結(jié)果。然后，將其導(dǎo)入到CFD，以上一次的流場分布為初始值，開始進行CFD仿真計算，實現(xiàn)CFD仿真數(shù)據(jù)的充分利用。

5．2 仿真并行化

雖然多級推力發(fā)動機的噴管在實際工作中的各個工作段是依次進行的，但由于將各個工作段的流場看成是定常流，因此各個工作段的流場是相對獨立的，對其進行的CFD仿真同樣是相互解偶的。因此，可采用并行化的方式，同時對不同工作段的流場進行CFD仿真。在生成網(wǎng)格后，分別導(dǎo)入不同工作段CFD仿真程序，按照不同的工作段的工況，分別設(shè)置不同的入口條件，之后就可同時進行仿真計算。這樣充分利用了多核或多機的計算資源，提高了計算效率。

5．3 優(yōu)化計算流程

按照圖4的流程圖，在Matlab平臺上完成優(yōu)化模型、流程控制和策略的實現(xiàn)，并通過成熟的商業(yè)計算流體力學軟件Fluent平臺搭建CFD模型，實現(xiàn)了最終的自動化的優(yōu)化計算流程。

圖4 優(yōu)化計算流程圖Fig．4 Flow chart of optimization calculation

6 結(jié)果和討論

6．1 問題描述

以某空空導(dǎo)彈固體火箭發(fā)動機長尾噴管的優(yōu)化問題為例，來驗證方法的有效性。

該發(fā)動機是一兩級推力的固體火箭發(fā)動機，長尾噴管的前長尾部分由于結(jié)構(gòu)限制無法變動，需要對尾部的收縮擴張噴管進行型面優(yōu)化，進一步挖掘工作潛力，提高噴管性能。

發(fā)動機噴管入口壓力為兩級壓力，分別為8．798 MPa和 4．719 MPa。要求噴管長度保持不變，并受到發(fā)動機及導(dǎo)彈結(jié)構(gòu)外形的嚴格約束，以此為約束條件，進行優(yōu)化計算。工作高度為地面高度，環(huán)境壓力為 0．101 325 MPa。

由于噴管的長度為固定長度，因此噴管的控制變量變?yōu)?個，即

6．2 計算結(jié)果

在8核Xeon計算服務(wù)器上，經(jīng)過64次的優(yōu)化迭代步，約54個機時完成了優(yōu)化計算，得到了優(yōu)化型面，如圖5所示。與原型面和優(yōu)化型面對應(yīng)的性能預(yù)測結(jié)果和實際試驗結(jié)果如表2所示。

其中，原型面噴管擴張段的二次曲線方程(原點為軸線上距出口138．35 mm處)為

對比優(yōu)化型面噴管擴張段的二次曲線方程(原點為軸線上距出口140．4 mm處)為

圖5 優(yōu)化型面和原型面的對比圖Fig．5 Comparison of optimized and original profile

表2 噴管型面性能優(yōu)化結(jié)果Table 2 Optimization result of performance for the nozzle profiles

根據(jù)表2的結(jié)果對比，對比原型面實際總沖預(yù)測值和試驗實測值，誤差只有0．18%，說明文中的CFD數(shù)值仿真模型的性能預(yù)測精度很高，噴管性能預(yù)測方法和結(jié)果可信。相比于CFD模型誤差，噴管的總沖優(yōu)化幅度達到了1．63%，其幅度遠比仿真誤差大，可認為優(yōu)化是基本可靠的，優(yōu)化結(jié)果是可信的。

總覽圖5和表2結(jié)果可知，原型面和優(yōu)化型面從兩方面提高了噴管的性能:

(1)從理想比沖和理想總沖的對比數(shù)據(jù)可知，通過調(diào)整噴管的擴張比，使得一級理想比沖降低，二級理想比沖增加，雖然二級比沖增加的幅度小于一級比沖，但由于二級工作段所消耗的推進劑質(zhì)量更多，使得理想總沖有0．08%小幅提高，這說明了優(yōu)化結(jié)果更好地平衡了噴管兩個工作階段的膨脹損失，減少了總的膨脹損失。

(2)從噴管效率的對比數(shù)據(jù)可知，通過優(yōu)化型面曲線，減少了噴管的摩擦損失等非膨脹損失，使噴管性能提高了1．55%，是優(yōu)化型面的性能提高的主要原因。

可看出，有效提高了噴管的總沖性能，達到1．63%，這對空空導(dǎo)彈的性能提升來說是十分重要的。

對比最終流場計算結(jié)果如表3和圖6所示。從中可看到，一級和二級工況的邊界總壓(最低總壓)都有顯著提高，而總壓的減少主要是由摩擦引起的，這說明優(yōu)化型面相對原型面顯著減少了噴管的摩擦損失。

通過對比兩相平衡流和兩相流的計算結(jié)果，可計算得到噴管原型面與優(yōu)化型面的兩相流損失，如表4所示。根據(jù)表2的噴管效率結(jié)果分別可計算出原型面和優(yōu)化型面噴管總損失為5．83%和4．37%。對比表4的結(jié)果可知，優(yōu)化型面相對原型面的兩相流損失雖然有0．14%的增加，從顆粒相體積分數(shù)分布圖可看到，如圖7所示，優(yōu)化型面沿壁面的顆粒分布更為集中，但相對噴管總損失的降低幅度1．55%，其影響很小。在該噴管的效率損失因素中，兩相流損失對噴管型面改變的敏感性相對較低。這說明本文所采用的優(yōu)化過程中，先后基于兩相平衡流和兩相流性能估算的分步優(yōu)化策略是有效的，并在保證結(jié)果精度的同時，降低了計算復(fù)雜度，提高了計算效率。

圖6 原型面和優(yōu)化型面的靜壓分布Fig．6 Static pressure distributions of original and optimized profile

表3 噴管內(nèi)最高與最低氣相總壓比較Table 3 Comparison of highest and lowest total pressure of gas phase in nozzle

表4 兩相流損失對比Table 4 Comparison of two-phase loss %

圖7 原型面和優(yōu)化型面的顆粒相體積分數(shù)分布Fig．7 Particle phase volume fraction distributions of original and optimized profile

7 結(jié)論

針對多級推力固體火箭發(fā)動機的特點，建立了可靠的CFD仿真模型，對不同工作段的噴管內(nèi)流場進行仿真，通過綜合各工作段噴管內(nèi)流場參數(shù)，計算得到發(fā)動機的實際總沖和噴管工作效率，實現(xiàn)了一種精確預(yù)估噴管性能的方法。

通過合理應(yīng)用CFD模型、動態(tài)網(wǎng)格、并行化仿真等技術(shù)和策略，將基于CFD仿真預(yù)估噴管性能的方法用于噴管型面的優(yōu)化設(shè)計中，解決了計算規(guī)模過大、計算不穩(wěn)定等問題，能高效穩(wěn)定、可靠地對多級推力的固體火箭發(fā)動機軸對稱噴管型面進行優(yōu)化計算。為固體火箭發(fā)動機設(shè)計人員深入地挖掘噴管潛能和提高噴管的設(shè)計水平，提供了切實可行和有效的方法。

考慮到基于數(shù)值仿真的優(yōu)化設(shè)計方法的通用性，本文對探索其他部件或領(lǐng)域的基于數(shù)值仿真的優(yōu)化設(shè)計計算方法也具有參考價值。

［1］ 方丁酉．兩相流噴管擴散段型面優(yōu)化計算［J］．航空動力學報，1986，1(2):137-140．

［2］ Cai G，F(xiàn)ang J，Xu X，et al．Performance prediction and optimization for liquid rocket engine nozzle［J］．Aerospace Science and Technology，2007，11(2-3):155-162．

［3］ Wang X，Damodaran M．Optimal three-dimensional nozzle shape design using CFD and parallel simulated annealing［J］．Journal of Propulsion and Power(0748-4658)，2002，18(1):217-221．

［4］ 李宜敏．固體火箭發(fā)動機原理［M］．北京:國防工業(yè)出版社，1985:308．

［5］ 劉海峰．單室雙推力固體火箭發(fā)動機裝藥［J］．上海航天，1990(3):5-9．

［6］ ANSYS Inc．Ansys fluent theory guide［Z］．ANSYS Inc．，2009．

［7］ Lun C K K，Savage S B，Jeffrey D J，et al．Kinetic theories for granular flow:inelastic particles in couette flow and slightly inelastic particles in a general flowfield［J］．Journal of Fluid Mechanics，1984，140:223-256．

［8］ Wen C Y，Yu Y H．Mechanics of fluidization［J］．Chemical Engineering Progress Symposium Series，1966，62(67):100-111．

［9］ Ghobadian A，Vasquez S A．A general purpose implicit coupled algorithm for the solution of eulerian multiphase transport equation:international conference on multiphase flow［Z］．Leipzig，Germany:University of Florida，2007．

［10］ 楊丹．固體火箭發(fā)動機氣-固兩相流的數(shù)值模擬［D］．哈爾濱工程大學，2006．

［11］ Waltz R A，Morales J L，Nocedal J，et al．An interior algorithm for nonlinear optimization that combines line search and trust region steps［J］．Mathematical Programming，2006，107(3)．

［12］ Nocedal J，Wright S J．Numerical optimization［M］．2ed．Springer Verlag，2006．

［13］ Thompson J F，Warsi Z U A，Mastin C W．Numerical grid generation:foundations and applications［M］．New York:North-Holland:Elsevier Science Pub．Co．，1985．