基于支持向量回歸模型的彈用沖壓發動機性能預測及優化

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

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

0 引言

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

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

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

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

1 發動機幾何模型與網格劃分

1.1 幾何模型

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

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

1.2 網格劃分

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

圖2 沖壓發動機網格Fig.2 Meshing of simulation cases

2 沖壓發動機內彈道計算模型

2.1 基本假設

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

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

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

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

2.2 計算模型

2.2.1 控制方程與湍流模型

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

(1)

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

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

(2)

(3)

式中:d為流動維數;Gk、Gω分別為湍動能和比耗散率的速度梯度;Yk和Yω分別為關于k和ω的湍流耗散項;Dω為交叉擴散項;Гk和Гω分別為關于k和ω的有效擴散系數。Transition SST模型的間歇因子和當地邊界層動量厚度雷諾數輸運方程見文獻[18]。

2.2.2 輻射模型

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

2.2.3 燃燒模型

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

C4H6+5.5O2→4CO2+3H2O。

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

(4)

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

(5)

表2 C4H6主要物性參數Table 2 Main parameters of C4H6

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

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

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

(6)

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

(7)

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

(8)

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

2.3 網格收斂性分析

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

表3 網格收斂性分析Table 3 Grid convergence analysis

將前5個網格與108萬網格的相比,發動機推力偏差分別為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%。在保證精度的情況下,為了盡可能地提高計算效率,選擇網格數目為72萬。

3 沖壓發動機工作性能優化模型

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

圖3 沖壓發動機優化流程Fig.3 Optimization workflow of theramjet

具體操作流程如下:

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

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

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

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

3.1 SVR模型

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

(9)

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

(10)

(11)

(12)

相應的預測函數變為

(13)

模型的誤差可表示為

(14)

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

3.2 多目標優化模型

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

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

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

(15)

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

3.3 帶有罰函數的NSGA-Ⅱ

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

(16)

(17)

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

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

(18)

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

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

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

4 數值模擬結果與分析

4.1 計算模型可靠性檢驗

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

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

圖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 優化結果分析

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

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

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

M=kt+λkI

(19)

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

圖9 λ=1時優化后的沖壓發動機結構Fig.9 Optimizedramjet shape when λ=1

表5 優化后發動機結構參數

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

圖10 發動機中C4H6的質量分數Fig.10 C4H6 mass fraction of the ramjet

圖11 發動機中O2的質量分數Fig.11 O2 mass fraction of the ramjet

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

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

圖13 優化后沖壓發動機的溫度和流線圖Fig.13 Temperature contour and streamline of the optimized ramjet

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

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

5 結論

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

表6 優化前后發動機性能指標Table 6 Performance of the ramjet

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

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

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