張寧, 史金光, 王中原, 趙新新

(南京理工大學能源與動力工程學院, 江蘇 南京 210094)

0 引言

將炮彈與固體燃料沖壓發(fā)動機相結合,可有效提高其射程,最大增程率可達70%,因此沖壓增程炮彈受到諸多國家的重視。對沖壓增程炮彈的研究始于20世紀70年代初期,此后美國、瑞典、南非、以色列等國家在沖壓增程理論及實驗方面取得了一些進展,提出了多種彈形結構,并就其空氣動力學及發(fā)動機燃燒特性開展了大量研究。文獻[1]首先對75 mm旋轉沖壓增程穩(wěn)定彈進行了研究,該彈采用皮托式進氣道,初速馬赫數(shù)為4.3,射程為12 km,實現(xiàn)了自點火和穩(wěn)定飛行,其后又研制出了采用中心體進氣道的203 mm尾翼穩(wěn)定彈,射程為60 km。文獻[2]進行了40 mm沖壓助推防空炮彈的研究,該彈的速度馬赫數(shù)約為4.3,燃燒時間約為2～3 s。文獻[3]開展了 155 mm 沖壓增程炮彈的實彈射擊,初速為 900 m/s,射程約為55 km。文獻[4]成功完成了使用沖壓發(fā)動機維持炮彈速度的飛行實驗。但這些研究大多都建立在旋轉彈的基礎上,提高射程的同時將會降低射擊精度[5],且也難以調整炮彈的飛行速度、高度與姿態(tài),不能控制其飛行彈道。為解決上述問題,文獻[6]采用簡易控制技術研制了155 mm沖壓增程制導炮彈,該炮彈頭部有控制翼,尾部有可展開尾翼,最大巡航速度馬赫數(shù)約為3,射程可達100 km,但加裝控制機構會進一步壓縮沖壓炮彈的內部空間,減小沖壓發(fā)動機體積,這將使發(fā)動機內燃氣速度增加、停留時間縮短,導致發(fā)動機推力、比沖和燃燒效率下降等問題[7]。

因此,有必要在考慮彈體幾何約束的情況下,優(yōu)化沖壓發(fā)動機結構,提高其工作性能。目前使用的優(yōu)化方法主要有兩種,一種是基于梯度的方法,另一種是啟發(fā)式算法。為了處理多目標優(yōu)化問題,前者需要通過一系列目標權重組合來定義全局目標函數(shù)[8],這種方法對假定的權重系數(shù)非常敏感,如果選擇了不合適的權重,可能會丟失一些最優(yōu)解。對于后者,目前使用最廣泛的方法是帶精英策略的非支配排序遺傳算法(NSGA-Ⅱ)[9],其具有計算速度快、解集收斂性好和能在單次優(yōu)化中生成帕累托前沿的優(yōu)點。

但在優(yōu)化過程中,計算流體力學(CFD)求解器通常必須被大量調用才能獲得最優(yōu)的解決方案,帶來優(yōu)化難度大、時間長等問題。于是,代理模型方法應運而生,并逐步受到重視,其基本思想是用一個簡單的逼近函數(shù)近似替代高精度求解器。較常用的代理模型有響應面[10]、徑向基函數(shù)[11]、人工神經網(wǎng)絡[12]和Kriging[13]等方法,雖然它們都具有較好的預測能力,但為了獲得較高的精度和泛化能力,必須使用大量的訓練樣本,從而削弱了代理模型的優(yōu)勢。支持向量回歸(SVR)模型[14]是基于結構最小化原則的一種機器學習方法,其核心思想是,基于Mercer核展開定理,將樣本空間映射到Hilbert空間中,并在其中應用線性方法來解決非線性回歸問題,保證了良好的泛化能力。因此,其在解決小樣本、非線性及高維問題中表現(xiàn)出了諸多優(yōu)勢[15]。

同時,由于沖壓增程制導炮彈體積較小,為了保證炮彈的威力、穩(wěn)定性和操縱性,通常僅能給發(fā)動機提供有限的安裝空間,限制了其長度、內徑和最大進氣面積等結構參數(shù)。因此,本文提出了一種與彈內有限空間適配的發(fā)動機性能預測及優(yōu)化方法,即基于帶轉捩的剪切應力輸運(TransitionSST)和渦概念耗散(EDC)方程,建立了內彈道計算模型;而后,基于SVR方法構建了性能預測模型,并結合NSGA-Ⅱ對發(fā)動機結構進行了優(yōu)化。所得結果可為彈用沖壓發(fā)動機的性能預測及結構優(yōu)化設計方法提供參考。

1 發(fā)動機幾何模型與網(wǎng)格劃分

1.1 幾何模型

本文以某沖壓增程彈為研究對象[5],圖1為其常用的發(fā)動機結構[16-17],總長為600 mm,入口直徑為35 mm;燃燒室、補燃室長度分別為247 mm和270 mm,裝藥內徑為75 mm;噴管喉部直徑為36 mm。

圖1 彈用沖壓發(fā)動機幾何模型Fig.1 Structure of SFRJ for projectile

1.2 網(wǎng)格劃分

圖2為發(fā)動機的網(wǎng)格劃分情況,在近壁面處加密以保證其附近參數(shù)的準確性。為模擬炮彈在海平面以馬赫數(shù)2.5飛行的工況,發(fā)動機入口的空氣質量流率為1.4 kg/s、溫度為540 K,總壓為1.0 MPa;噴管出口為壓力出口;壁面絕熱。

圖2 沖壓發(fā)動機網(wǎng)格Fig.2 Meshing of simulation cases

2 沖壓發(fā)動機內彈道計算模型

2.1 基本假設

為簡化仿真過程,獲得發(fā)動機的主要工作性能,做出如下假設:

1) 燃氣可近似視為理想氣體;

2) 燃料內壁為氣固耦合交界面,外壁為絕熱壁;

3) 端羥基聚丁二烯(HTPB)推進劑的熱解產物為1.3-丁二烯單質(C4H6)。

2.2 計算模型

2.2.1 控制方程與湍流模型

帶化學反應的軸對稱雷諾時均Navier-Stokes(N-S)方程如下:

(1)

式中:x、r分別為軸向和徑向坐標;ρ為密度;p為壓強;u、v分別為軸向和徑向速度;μl、μt分別為分子和湍流黏性系數(shù);Prl和Prt分別為分子和湍流普朗特數(shù);Cp為定壓熱容;Cv為定容熱容;T為溫度;mj為組分質量分數(shù);Rj、Dj分別為化學組分j的反應和擴散速率。

空間離散采用2階迎風型矢通量分裂格式。由于發(fā)動機的內流場包含由尖銳幾何形狀引起的轉捩區(qū)域[7],湍流模型選用Transition SST模型,其湍動能k以及比耗散率ω的輸運方程為

(2)

(3)

式中:d為流動維數(shù);Gk、Gω分別為湍動能和比耗散率的速度梯度;Yk和Yω分別為關于k和ω的湍流耗散項;Dω為交叉擴散項;Гk和Гω分別為關于k和ω的有效擴散系數(shù)。Transition SST模型的間歇因子和當?shù)剡吔鐚觿恿亢穸壤字Z數(shù)輸運方程見文獻[18]。

2.2.2 輻射模型

本文選用離散坐標輻射模型[19]模擬發(fā)動機中的輻射傳熱。

2.2.3 燃燒模型

空氣進入燃燒室后與C4H6發(fā)生如下化學反應[20]:

C4H6+5.5O2→4CO2+3H2O。

化學反應速率采用渦耗散模型計算。C4H6由藥柱表面熱解產生,因此在計算過程中,需在燃料表面與流體域交界面的第1層網(wǎng)格上加質。HTPB的熱解速率與燃料表面溫度有關,服從Arrhenius公式:

(4)

式中:A為指前因子;Ea為活化能;R為氣體常數(shù);Tw為燃料內壁溫度,可由氣相-固相分界面上的能量平衡方程求解[21]:

(5)

表2 C4H6主要物性參數(shù)Table 2 Main parameters of C4H6

根據(jù)文獻[16-17],燃料的燃速(通常不超過1 mm/s)相較于發(fā)動機內氣流速度(一般為100 m/s)較慢,因此,可忽略燃面退移對流場的影響。

2.2.4 推力與比沖、燃燒效率模型

由動量定理,發(fā)動機推力F為出口截面氣流與迎面氣流的沖量差,即

(6)

比沖可定義為單位質量燃料所產生的推力,即

(7)

根據(jù)文獻[21],燃燒效率可定義為

(8)

式中:YCO2、YC4H6表示混合氣體中各成分的質量分數(shù);MCO2、MC4H5表示混合氣體中各成分的摩爾質量;S為截面面積。

2.3 網(wǎng)格收斂性分析

由于在固定幾何形狀下,發(fā)動機的比沖僅與其推力與燃速有關,可只選取初始沖壓發(fā)動機的推力、燃燒效率和平均燃速作為網(wǎng)格獨立性檢驗的指標。選取網(wǎng)格的總數(shù)約為14、21、32、48和72萬個,收斂結果如表3所示。

表3 網(wǎng)格收斂性分析Table 3 Grid convergence analysis

將前5個網(wǎng)格與108萬網(wǎng)格的相比,發(fā)動機推力偏差分別為38.82%、20.67%、7.47%、2.78%、0.28%,燃燒效率偏差為9.16%、7.33%、3.32%、2.31%、0.97%,平均燃速偏差為14.51%、12.16%、3.99%、2.53%、0.73%。在保證精度的情況下,為了盡可能地提高計算效率,選擇網(wǎng)格數(shù)目為72萬。

3 沖壓發(fā)動機工作性能優(yōu)化模型

受制導炮彈結構限制,其所用的沖壓發(fā)動機通常較小,同時發(fā)動機結構對其性能有極大影響[22]。例如,當發(fā)動機總長一定時燃燒室越長,則推進劑裝填量越大,能產生更大的推力,但比沖較低;反之,發(fā)動機的比沖較大而推力較小[23]。因此,有必要在考慮炮彈幾何約束的情況下,優(yōu)化發(fā)動機結構以提高其推進性能。圖3為本文建立的優(yōu)化工作流程,包含預處理、內彈道計算和多目標優(yōu)化3個模塊。

圖3 沖壓發(fā)動機優(yōu)化流程Fig.3 Optimization workflow of theramjet

具體操作流程如下:

1) 使用拉丁超立方抽樣(LHS)方法[24]在設計空間中生成訓練和測試數(shù)據(jù)集。

2) 將訓練和測試集中的每個樣本代入內彈道模型進行計算,得到發(fā)動機的推力、比沖、燃燒效率等指標。

3) 建立SVR模型并通過測試集進行精度校驗,若不滿足,則對其參數(shù)進行尋優(yōu)更新。

4) 基于多目標優(yōu)化模型得到最終的優(yōu)化解。

3.1 SVR模型

f(x)=w·φ(x)+b

(9)

式中:w、b分別為權向量和閾值;φ(x)為設計變量的組合函數(shù),w·φ(x)為w與φ(x)的內積,并且滿足結構風險最小化原理。在精度ε下,f(x)能夠估計每組設計變量對應的響應值,即為所需的近似函數(shù),其參數(shù)求解可以轉化成求解凸優(yōu)化問題:

(10)

(11)

(12)

相應的預測函數(shù)變?yōu)?/p>

(13)

模型的誤差可表示為

(14)

式中:Nc為測試集中的樣本數(shù)目;為SVR估計響應值。若模型精度不滿足要求,則可使用遺傳算法對參數(shù)C、γ尋優(yōu)更新。

3.2 多目標優(yōu)化模型

圖4為沖壓發(fā)動機的幾何模型,其中,L為沖壓發(fā)動機總長,din為進氣道入口截面的半徑;hc、lc分別為燃燒室突擴臺階高度和長度;hp為隔板高度;lm、dm分別為補燃室的長度和內徑;dr為尾噴管的喉部半徑,ls為收斂段長度,le為擴張段長度,dout為尾噴管的出口半徑,本文選用hc、lc、hp、lm、dm、dr、ls、le、dout這9個參數(shù)作為優(yōu)化設計變量。

圖4 沖壓發(fā)動機幾何模型Fig.4 Geometric model of the ramjet

在發(fā)動機性能優(yōu)化過程中,除了要保證燃料的混合效率和燃燒效率外,還需要考慮彈上其余裝置對空間的需求以及與進氣道的適配性。例如,由于發(fā)動機的補燃室上方帶有戰(zhàn)斗部,為保證彈藥威力,燃燒室和噴管的長度、補燃室內徑都不宜過大。根據(jù)以上分析,沖壓發(fā)動機的工作性能優(yōu)化模型建立如下:

(15)

式中:kt、kI分別為推力、比沖的歸一化比值;F0和I0分別為具有初始沖壓發(fā)動機的推力和比沖值;Lmax為發(fā)動機最大長度;hmax為突擴臺階最大高度;dmax為補燃室最大內徑。式(15)約束條件中有:設計變量的取值范圍,發(fā)動機尺寸與炮彈可用空間的適用性;沖壓發(fā)動機性能約束,以確保大部分燃料能量可以在發(fā)動機中釋放;適配性約束,確保發(fā)動機產生的壓力不會超過氣道所能抵抗的最大背壓pmax。

3.3 帶有罰函數(shù)的NSGA-Ⅱ

由于沖壓發(fā)動機工作性能優(yōu)化問題是帶有約束的多目標優(yōu)化問題,可采用NSGA-Ⅱ[9]與罰函數(shù)相結合的方法進行求解,即將原適應度函數(shù)與Static Hoffmeister(SH)懲罰項組合,構造增廣適應度函數(shù):

(16)

(17)

式中:f0(x)為原適應度函數(shù);?為罰因子;gz(x)為第z個約束函數(shù);H(gi(x))在滿足約束條件時取0,否則取1。為了提高計算效率,根據(jù)范數(shù)相容性原理,可將式(17)改寫為

fitness(x)=f0(x)+?H(gmax)|gmax|

(18)

式中:gmax為gi(x)的最大值。

帶有罰函數(shù)的NSGA-Ⅱ算法流程如圖5所示。

圖5 帶有罰函數(shù)的NSGA-Ⅱ流程Fig.5 NSGA-Ⅱ flow witha penalty function

4 數(shù)值模擬結果與分析

4.1 計算模型可靠性檢驗

4.1.1 內彈道計算模型可靠性驗證

依照文獻[25]所述方法,本文對所用的內彈道計算模型進行了可靠性驗證。針對冷流實驗[26],獲得的燃燒室中心線軸向速度va與實驗結果基本一致,如圖6所示,而回流區(qū)長度較實驗值約差0.98%,表明該計算模型可較好地模擬發(fā)動機內氣流的流動過程。

圖6 燃燒室中心軸向速度Fig.6 Mean axial velocity alongthe combustor centerline

表4 計算結果與實驗結果的對比Table 4 Comparison between the calculated results and experimental results

4.1.2 SVR模型可靠性驗證

為檢驗代理模型的預測精度,采用100個樣本點建立了預測模型,并使用另外40個樣本點進行計算,與高精度求解器結果間的誤差對比如圖7所示。圖7中,δt表示推力誤差,δI表示比沖誤差,模型預測結果的平均誤差為1.24%,最大誤差為2.71%,不超過3%,可靠性較好。

圖7 預測值與計算值對比Fig.7 Comparison between predicted and calculated values

4.2 優(yōu)化結果分析

圖8顯示了NSGA-Ⅱ獲得的非支配結果,可以看到較清晰的帕累托前沿分布,比沖隨著推力的增加而減小,表明在此工況下沖壓發(fā)動機的優(yōu)化設計必須在兩個目標函數(shù)間進行折衷。

圖8 沖壓發(fā)動機多目標優(yōu)化設計的帕累托解集Fig.8 Pareto solution set for multi-objective optimization design

因此,本文設計了一個決策函數(shù),在解集中選擇一個合適的設計點,即

M=kt+λkI

(19)

式中:λ表示權重系數(shù),可以根據(jù)目標的重要性進行調整,在此分析中該值取為1?；诖?獲得的最佳沖壓發(fā)動機形狀如圖9所示,參數(shù)如表5所示。由圖9可以發(fā)現(xiàn),優(yōu)化后的沖壓發(fā)動機燃燒室長度減小,補燃室長度增加、內徑減小,一方面可以增加燃氣的摻混效果,提高發(fā)動機的燃燒效率;另一方面則可使補燃室上方戰(zhàn)斗部的裝藥量增加,保證了彈藥威力。

圖9 λ=1時優(yōu)化后的沖壓發(fā)動機結構Fig.9 Optimizedramjet shape when λ=1

表5 優(yōu)化后發(fā)動機結構參數(shù)

圖10、圖11出了優(yōu)化前后沖壓發(fā)動機中C4H6和O2的分布情況,可以發(fā)現(xiàn),二者組分分布情況一致,C4H6都主要集中于發(fā)動機壁面附近,而氧氣則主要位于發(fā)動機的中心區(qū)域,火焰層在富氧區(qū)與燃料內壁間形成,因此燃氣的混合效果較差,降低了發(fā)動機的燃燒速率。因此,部分燃料將在補燃室中繼續(xù)和氧氣反應。但在優(yōu)化后發(fā)動機尤其是補燃室內C4H6剩余較少。這是因為其補燃室更長,提高了燃料的停留時間,使其與空氣反應更為充分;同時燃燒室內的湍流動能隨著突擴臺階高度的增加而增大[17],提高了室內燃料與空氣的摻混和燃燒效率。

圖10 發(fā)動機中C4H6的質量分數(shù)Fig.10 C4H6 mass fraction of the ramjet

圖11 發(fā)動機中O2的質量分數(shù)Fig.11 O2 mass fraction of the ramjet

圖12 初始沖壓發(fā)動機的溫度和流線圖Fig.12 Temperaturecontour and streamline of the initial ramjet

圖12、圖13給出了優(yōu)化前后沖壓發(fā)動機的溫度和流線圖。分析流線圖可知,二者流動情況大致相同,但優(yōu)化后,燃燒室內回流區(qū)長度增加、補燃室內回流區(qū)長度減小,這是因為其突擴臺階高度增加,補燃室內徑減小,改變了原有結構的后臺階高度,從而改變了其后部漩渦的再附長度。由二者的溫度圖像可知,發(fā)動機優(yōu)化后,其燃燒室后部和補燃室內的溫度更高,這是因為其突擴臺階高度的增加,提高了燃燒室內的湍流動能,使得燃燒室內燃氣摻混和燃燒效率增加,同時,其較長的補燃室保證了燃氣在其中充分的燃燒,這也與圖10、圖11的C4H6和O2的分布情況相符。

圖13 優(yōu)化后沖壓發(fā)動機的溫度和流線圖Fig.13 Temperature contour and streamline of the optimized ramjet

圖14 燃速沿軸向分布情況Fig.14 Regression rate along grain surface

表6比較了優(yōu)化前后的推力和比沖。由表6可知,優(yōu)化后發(fā)動機推力提高了24.22%。由圖14可知,初始沖壓發(fā)動機的燃速約為0.774 mm/s,而優(yōu)化沖壓發(fā)動機的燃速約為0.822 mm/s。后者雖然燃速較高,但由于燃料長度縮短了13.88%,其燃料質量流量與前者相近,同時優(yōu)化后燃燒效率約提高了12.02%,其比沖也提高了約20.28%。因此優(yōu)化后的沖壓式噴氣發(fā)動機具有更好的推進性能。

5 結論

本文提出了一種與彈內有限空間適配的沖壓發(fā)動機性能預測及優(yōu)化方法。首先,采用Transition SST和EDC方程,建立了其內彈道計算模型,并獲得了相應的流場結構與性能參數(shù);在此基礎上,基于SVR模型構建了發(fā)動機性能參數(shù)的預測模型,并與NSGA-Ⅱ結合對發(fā)動機結構進行了優(yōu)化。得到以下主要結論:

表6 優(yōu)化前后發(fā)動機性能指標Table 6 Performance of the ramjet

1) 在工況相同的情況下,本文方法計算所得的燃料平均燃速,與實驗結果相比,絕對誤差平均值不超過1.5%,表明所用模型精度較高,能較好地模擬發(fā)動機內的燃燒與流動過程。同時,本文構建代理模型的預測結果與高可信度模型的計算結果相差較小,最大相對誤差不超過3%,具有較高的精度,可以提升沖壓發(fā)動機性能計算與分析的效率。

2) 優(yōu)化前后流場結構與各組分分布情況大致相同,都在發(fā)動機燃燒室入口后臺階以及補燃室隔板后臺階處出現(xiàn)了回流區(qū),但優(yōu)化后的燃燒室內回流區(qū)更長;C4H6主要集中于燃燒室上部,而氧氣主要集中分布在燃燒室通道的中心區(qū)域,但優(yōu)化后的發(fā)動機尤其是補燃室內C4H6剩余較少。

3)發(fā)動機結構優(yōu)化后,燃燒室縮短了13.88%,補燃室長度增加了13.50%,空氣與燃料混合更充分,推力、比沖和燃燒效率分別增加了24.22%、20.28%、12.02%。利用該優(yōu)化設計方法設計的發(fā)動機性能更優(yōu),為有限空間內彈用沖壓發(fā)動機的設計提供了可行的方法與工具。