陳 偉，梁國(guó)柱

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

首先，從以下2個(gè)方面來(lái)簡(jiǎn)化CFD仿真模型，減少計(jì)算規(guī)模。

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

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

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

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

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

2 CFD仿真模型

2．1 基本假設(shè)

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

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

(2)燃?xì)庵械臍庀酁橥耆珰怏w，服從完全氣體狀態(tài)方程;

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

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

2．2 控制方程

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

連續(xù)方程為

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

動(dòng)量守恒方程為

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

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

阻力系數(shù)CD為

能量方程為

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

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

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

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

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

入口邊界條件為壓力入口，其壓力為各個(gè)階段壓力的平均值。對(duì)于兩相平衡流，入口總溫為燃燒室燃?xì)饪倻?。?duì)于氣固兩相流，氣相壓力和顆粒相均采用兩相平衡流的燃燒室燃?xì)饪倻刈鳛槿肟诳倻兀w粒相的入口體積分?jǐn)?shù)由推進(jìn)劑顆粒相質(zhì)量分?jǐn)?shù)和氣相與顆粒相的入口密度換算得到。

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

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

3 噴管性能預(yù)估

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

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

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

引入熱力計(jì)算得到的噴管入口的燃?xì)鈪?shù)，可直接計(jì)算得到各個(gè)工作段不計(jì)噴管損失的理想比沖，并進(jìn)一步計(jì)算得到理想總沖，即

由于噴管入口的燃?xì)鈪?shù)是給定的，因此計(jì)算噴管的總工作效率為實(shí)際比沖除以理想比沖，即

通過(guò)以上步驟計(jì)算得到噴管總工作效率η和實(shí)際總沖I指標(biāo)，就可對(duì)多級(jí)推力的固體火箭發(fā)動(dòng)機(jī)的噴管性能做出精確預(yù)估和評(píng)價(jià)。

4 優(yōu)化模型

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

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

4．2 數(shù)學(xué)模型

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

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

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

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

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

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

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

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

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

5 計(jì)算流程實(shí)現(xiàn)

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

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

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

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

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

5．2 仿真并行化

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

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

按照?qǐng)D4的流程圖，在Matlab平臺(tái)上完成優(yōu)化模型、流程控制和策略的實(shí)現(xiàn)，并通過(guò)成熟的商業(yè)計(jì)算流體力學(xué)軟件Fluent平臺(tái)搭建CFD模型，實(shí)現(xiàn)了最終的自動(dòng)化的優(yōu)化計(jì)算流程。

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

6 結(jié)果和討論

6．1 問(wèn)題描述

以某空空導(dǎo)彈固體火箭發(fā)動(dòng)機(jī)長(zhǎng)尾噴管的優(yōu)化問(wèn)題為例，來(lái)驗(yàn)證方法的有效性。

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

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

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

6．2 計(jì)算結(jié)果

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

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

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

圖5 優(yōu)化型面和原型面的對(duì)比圖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é)果對(duì)比，對(duì)比原型面實(shí)際總沖預(yù)測(cè)值和試驗(yàn)實(shí)測(cè)值，誤差只有0．18%，說(shuō)明文中的CFD數(shù)值仿真模型的性能預(yù)測(cè)精度很高，噴管性能預(yù)測(cè)方法和結(jié)果可信。相比于CFD模型誤差，噴管的總沖優(yōu)化幅度達(dá)到了1．63%，其幅度遠(yuǎn)比仿真誤差大，可認(rèn)為優(yōu)化是基本可靠的，優(yōu)化結(jié)果是可信的。

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

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

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

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

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

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

圖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 兩相流損失對(duì)比Table 4 Comparison of two-phase loss %

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

7 結(jié)論

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

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

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

［1］ 方丁酉．兩相流噴管擴(kuò)散段型面優(yōu)化計(jì)算［J］．航空動(dòng)力學(xué)報(bào)，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ā)動(dòng)機(jī)原理［M］．北京:國(guó)防工業(yè)出版社，1985:308．

［5］ 劉海峰．單室雙推力固體火箭發(fā)動(dòng)機(jī)裝藥［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ā)動(dòng)機(jī)氣-固兩相流的數(shù)值模擬［D］．哈爾濱工程大學(xué)，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．