趙子熹，姜毅，賈啟明，李玉龍，牛鈺森

(1.北京理工大學(xué) 宇航學(xué)院， 北京 100081； 2.火箭軍研究院 系統(tǒng)工程研究所， 北京 100085)

0 引言

彈射技術(shù)在航空航天特別是導(dǎo)彈發(fā)射領(lǐng)域得到廣泛應(yīng)用。除電磁彈射外，其余導(dǎo)彈彈射方式均依賴工作介質(zhì)做功，以滿足設(shè)計的彈射指標(biāo)要求。根據(jù)工作介質(zhì)的類型，導(dǎo)彈彈射方式可分為火藥燃?xì)馐?、燃?xì)? 蒸汽混合式、壓縮氣體式、液力式和超臨界流體式等。

現(xiàn)階段相對成熟的火藥燃?xì)馐綇椛溲b置存在熱防護(hù)要求較高、發(fā)射時隱蔽性較差等缺點；壓縮氣體式彈射裝置則在設(shè)備復(fù)雜度、彈射有效載荷等方面表現(xiàn)出一定的劣勢。超臨界流體特別是超臨界二氧化碳一方面具有較高的工作壓強(qiáng)，可適配較大區(qū)間的彈射載荷需求；另一方面工作溫度相對較低，在彈射過程中對導(dǎo)彈的影響較小，同時還具有比內(nèi)能相對較高、發(fā)射時光熱信號強(qiáng)度低、清潔環(huán)保等優(yōu)點，逐步進(jìn)入研究者的視野。

對基于超臨界流體的彈射裝置進(jìn)行研究，一個關(guān)鍵在于正確表征流體的熱物性參數(shù)。超臨界流體的壓強(qiáng)一般較高，經(jīng)典的理想氣體假設(shè)已不再適用，可按照真實氣體進(jìn)行處理。Ren等基于改進(jìn)的維里方程構(gòu)建以高壓空氣為工質(zhì)的彈射動力學(xué)方程，對比了基于真實和理想氣體假設(shè)的仿真結(jié)果，得到了真實氣體效應(yīng)對彈射動力學(xué)特性有顯著影響的結(jié)論。姚琳等基于Peng-Robinson方程構(gòu)建了兩級提拉式單側(cè)彈射裝置內(nèi)彈道模型并通過縮比實驗驗證了仿真結(jié)果；在建模仿真的基礎(chǔ)上以氣源體積為優(yōu)化目標(biāo)建立優(yōu)化模型，有效提升了發(fā)射裝置機(jī)動性。Yao等對基于超臨界二氧化碳的無人機(jī)彈射裝置進(jìn)行了研究，在構(gòu)建氣- 液兩相彈射動力學(xué)方程的基礎(chǔ)上探討了高壓室初始狀態(tài)等參數(shù)對彈射過程的影響，發(fā)現(xiàn)與空氣等傳統(tǒng)工質(zhì)相比，超臨界二氧化碳作為一種新型工作介質(zhì)具有較高的工作效能。

前述文獻(xiàn)在建立動力學(xué)模型時，均采用具有高壓室- 噴管/閥門- 低壓室形式的彈射裝置，在處理噴管流動時引入具有解析形式的流量表達(dá)式，并未討論真實氣體效應(yīng)對噴管流動造成的影響。由于超臨界流體流經(jīng)噴管時仍具有較高壓強(qiáng)，不能被簡單視為理想氣體，有必要作為真實氣體處理。同時部分文獻(xiàn)中的彈射內(nèi)彈道優(yōu)化模型主要采用的是單目標(biāo)優(yōu)化思想，難以從多個維度反映彈射裝置效能。

本文針對使用超臨界二氧化碳作為工質(zhì)的導(dǎo)彈彈射裝置，基于Soave-Redlich-Kwong(S-R-K)狀態(tài)方程建立了彈射內(nèi)彈道動力學(xué)模型，特別是考慮了工質(zhì)經(jīng)噴管由高壓室向低壓室流動過程的真實氣體效應(yīng)。在內(nèi)彈道模型的基礎(chǔ)上構(gòu)建多目標(biāo)優(yōu)化模型，對內(nèi)彈道參數(shù)進(jìn)行了優(yōu)化設(shè)計。計算和優(yōu)化結(jié)果對超臨界二氧化碳在彈射系統(tǒng)中的應(yīng)用具有一定的參考意義。

1 基于超臨界二氧化碳的彈射動力學(xué)模型

1.1 模型假設(shè)

導(dǎo)彈彈射裝置采用垂直發(fā)射裝置，主要由高壓室、噴管和低壓室三部分組成，結(jié)構(gòu)如圖1所示。在彈射過程開始前，超臨界二氧化碳存儲于高壓室中，噴管處于關(guān)閉狀態(tài)。彈射開始后，噴管打開，超臨界二氧化碳進(jìn)入低壓室，推動導(dǎo)彈運動。

圖1 彈射裝置結(jié)構(gòu)簡圖Fig.1 Schematic diagram of the catapult

在構(gòu)建內(nèi)彈道動力學(xué)模型前，給定如下基本假設(shè)：

1)彈射過程視為單相過程，即將彈射裝置任一組成部分中的工質(zhì)作為均一相處理，不考慮多相平衡共存的情況。

2)根據(jù)文獻(xiàn)[5]的實驗及數(shù)值仿真結(jié)果，噴管內(nèi)超臨界二氧化碳流動存在非平衡凝結(jié)等現(xiàn)象，增加了建模和計算的復(fù)雜性。為簡化問題，引入準(zhǔn)平衡思想，認(rèn)為超臨界二氧化碳在噴管中的流動由一系列平衡狀態(tài)構(gòu)成，忽略非平衡過程帶來的影響。在假設(shè)1的基礎(chǔ)上認(rèn)為噴管內(nèi)工質(zhì)流動屬于單相、平衡過程。

3)工質(zhì)在噴管內(nèi)流動視為一維等熵絕熱流動，不考慮管內(nèi)激波、摩擦等對流動造成的影響。

4)工質(zhì)在噴管內(nèi)流動始終作為臨界流動處理，即噴管喉部流動馬赫數(shù)始終為1。

5)除噴管外高壓室和低壓室內(nèi)工質(zhì)狀態(tài)量在空間均勻分布，只隨時間變化。

6)導(dǎo)彈在低壓室的運動是無摩擦的，即忽略導(dǎo)彈適配器等因素的影響。

7)忽略系統(tǒng)與外界環(huán)境之間的能量交換。

1.2 真實氣體狀態(tài)方程

1.2.1 Soave-Redlich-Kwong狀態(tài)方程

理論建模中一般通過狀態(tài)方程引入真實氣體效應(yīng)。S-R-K狀態(tài)方程是Soave基于Redlich-Kwong方程提出的一種立方型狀態(tài)方程，以在多個應(yīng)用領(lǐng)域中體現(xiàn)出的精度而著稱。S-R-K方程可表示為

(1)

式中：為工質(zhì)壓強(qiáng)；為特定氣體常數(shù)；為工質(zhì)溫度；為工質(zhì)比體積；為物質(zhì)偏心因子；、為特定物質(zhì)常數(shù)，與特定物質(zhì)的關(guān)系為

、分別為物質(zhì)的臨界溫度和臨界壓強(qiáng)；(，)為與工質(zhì)溫度相關(guān)的修正函數(shù)，

(2)

122 超臨界二氧化碳熱物性參數(shù)建模

基于選定的狀態(tài)方程，構(gòu)建超臨界二氧化碳定壓比熱容、比焓等狀態(tài)量計算表達(dá)式。計算公式為

(3)

式中：為工質(zhì)定容比熱容；、和間的偏導(dǎo)數(shù)可根據(jù)狀態(tài)方程(1)式及偏導(dǎo)數(shù)關(guān)系得到

(4)

(5)

(6)

對于真實氣體或超臨界流體而言，不僅與溫度有關(guān)，還受到壓強(qiáng)的影響。根據(jù)Maxwell關(guān)系式及狀態(tài)方程表達(dá)式，有

(7)

設(shè)想工質(zhì)密度由0 kg/m變化至特定值的等溫過程，則可表示為

(8)

(9)

(10)

(=1,2，…，5)為系數(shù)，=2535 6，=8067 1×10，=-5706 1×10，=1793 9×10，=-1528 2×10。

(10)式的適用范圍為100～1 477 K。將(4)式～(10)式代入(3)式，即可計算特定壓強(qiáng)和溫度下工質(zhì)定壓比熱容。

比焓可表示為

(11)

式中：為工質(zhì)比內(nèi)能，對于真實氣體，為和的函數(shù)，

(12)

設(shè)想一個工質(zhì)密度由0 kg/m變化至的等溫過程，則的變化為

(13)

綜合(9)式、(10)式、(12)式、(13)式及(4)式，可得的表達(dá)式為

(14)

式中：為常數(shù)系數(shù)，=58817 9。結(jié)合(11)式、(14)式可得特定壓力及溫度下的工質(zhì)比焓計算式。

根據(jù)相同建模思想可獲得工質(zhì)比熵的計算式為

(15)

式中：為常數(shù)系數(shù)，=9895×10。

123 S-R-K狀態(tài)方程驗證

為確認(rèn)S-R-K狀態(tài)方程的精度，將基于該狀態(tài)方程計算的二氧化碳壓強(qiáng)、定壓比熱容、比焓和比熵等狀態(tài)量與美國國家標(biāo)準(zhǔn)與技術(shù)研究院(NIST)數(shù)據(jù)庫中的相關(guān)數(shù)據(jù)進(jìn)行對比，結(jié)果如圖2所示。使用的二氧化碳物性參數(shù)如表1所示。

表1 二氧化碳熱物性參數(shù)

以NIST為基準(zhǔn)時計算值平均相對誤差如表2所示。

表2 熱物性參數(shù)平均相對誤差

根據(jù)圖2及表2，在給定的溫度范圍內(nèi)當(dāng)二氧化碳密度小于500 kg/m時，基于S-R-K方程計算的壓強(qiáng)、定壓比熱容、比焓和比熵等物理量與NIST數(shù)據(jù)較為一致，平均相對誤差均維持在5%以內(nèi)。隨著密度增大，壓強(qiáng)和比熵計算誤差逐漸增大，接近或超過30%，表明S-R-K方程預(yù)測二氧化碳熱物性參數(shù)的能力有所減弱；定壓比熱容的計算誤差增加至約15%后減小，同時熱物性參數(shù)模型也相對準(zhǔn)確地預(yù)測了臨界點附近的參數(shù)躍變；比焓計算值的偏差始終較小，平均相對誤差始終保持在4%以內(nèi)。

圖2 基于S-R-K方程的二氧化碳熱物性參數(shù)計算值與NIST數(shù)據(jù)對比Fig.2 Comparison of the thermal physical results of carbon dioxide calculated by S-R-K equation with NIST data

綜合上述分析，認(rèn)為當(dāng)密度處于一定范圍內(nèi)(≤500 kg/m)時，S-R-K狀態(tài)方程預(yù)測的二氧化碳物性參數(shù)計算誤差相對較小，可以滿足后續(xù)內(nèi)彈道計算需求。

1.3 考慮真實氣體效應(yīng)的一維等熵噴管流動模型

1.3.1 一維等熵噴管流動模型

噴管作為連通高、低壓室的部件，其流動特性直接決定低壓室工質(zhì)狀態(tài)，進(jìn)而影響彈射效能。根據(jù)能量守恒原理，一維噴管流動的能量方程為

(16)

式中：為滯止?fàn)顟B(tài)比焓；為噴管中距入口距離處的工質(zhì)流動速度。以下標(biāo)t表示噴管喉部，則噴管流量可表示為

(17)

式中：為噴管喉部橫截面積；為滯止?fàn)顟B(tài)壓力；為滯止?fàn)顟B(tài)密度；為無量綱流動系數(shù)，

(18)

為噴管喉部密度，為噴管喉部的比焓。根據(jù)(17)式、(18)式即可確定噴管臨界流量。

132 一維等熵噴管流動的迭代解法

根據(jù)熵的定義式和Maxwell關(guān)系式，熵可表示為

(19)

基于等熵流動假設(shè)d=0，可得

(20)

真實氣體比焓可表示為溫度和壓強(qiáng)的函數(shù)：

(21)

將(20)式、(21)式改寫為如下迭代形式：

(22)

式中：下標(biāo)表示第個網(wǎng)格節(jié)點。根據(jù)噴管流動的定性分析，臨界流動下溫度沿流向降低，即d始終為負(fù)值。因此可以沿軸向?qū)姽芸臻g離散為若干網(wǎng)格節(jié)點，給定計算步長d，則溫度的迭代計算式為

=-1+d

(23)

迭代計算需要設(shè)置計算終止條件。彈射內(nèi)彈道問題對于噴管部分，主要關(guān)注噴管流量。根據(jù)臨界流動假設(shè)，流動達(dá)到臨界狀態(tài)后，流量不再隨噴管出入口壓力比的減小而增加。臨界流量在工程實踐中常用喉部流量表示。同時噴管內(nèi)馬赫數(shù)沿流向增加，在喉部始終為1。綜上，迭代終止準(zhǔn)則可設(shè)置為：當(dāng)網(wǎng)格點處馬赫數(shù)≥1時停止計算，最接近1的馬赫數(shù)對應(yīng)的網(wǎng)格節(jié)點可以表示噴管喉部。在確定喉部位置的基礎(chǔ)上，根據(jù)(17)式、(18)式計算的流量值即為噴管臨界流量。為計算馬赫數(shù)，首先根據(jù)聲速定義有

(24)

式中：表示聲速。等熵條件下壓強(qiáng)對比體積的偏導(dǎo)數(shù)根據(jù)狀態(tài)方程推導(dǎo)為

(25)

對于右端新增的等熵條件下溫度對比體積的偏導(dǎo)數(shù)，注意到熵還可表示為溫度和比體積的函數(shù)：

(26)

由等熵假設(shè)可得：

(27)

綜合(4)式、(25)式、(27)式并將其代入(24)式即得聲速計算式。在解出聲速的基礎(chǔ)上即可確定對應(yīng)網(wǎng)格點處馬赫數(shù)。

以上給出了一維等熵噴管流動迭代解法的基本思想。實際計算時，首先確定計算步長d，之后根據(jù)(22)式、(23)式進(jìn)行迭代求解。求得網(wǎng)格點處溫度、壓強(qiáng)和比焓后，密度由狀態(tài)方程(1)式確定，速度由能量方程(16)式確定。

133 噴管流動模型驗證

在構(gòu)建噴管流動模型的基礎(chǔ)上進(jìn)行網(wǎng)格無關(guān)性驗證，網(wǎng)格步長d取值范圍為-005～-0002 K。初始條件設(shè)置如表3所示，噴管喉部直徑為85 mm。

表3 網(wǎng)格無關(guān)性驗證初始條件

網(wǎng)格無關(guān)性驗證結(jié)果如圖3所示。

圖3 噴管流動模型網(wǎng)格無關(guān)性驗證Fig.3 Grid independence validation of the nozzle flow model

根據(jù)圖3可知，隨著計算步長絕對值的減小，流量和無量綱流動系數(shù)均呈現(xiàn)增長趨勢。若以計算步長-0.05 K作為對照基準(zhǔn)，則流量最大絕對差值為0.059 kg/s，對應(yīng)相對差值為0.029 8%；無量綱流動系數(shù)最大絕對差值為0.000 15，對應(yīng)相對差值為0.029 9%。驗證計算表明，計算步長對噴管流動流動模型計算結(jié)果影響較小，后續(xù)計算時可在綜合計算時間成本等因素基礎(chǔ)上在-0.03～-0.01 K的范圍內(nèi)設(shè)置步長。

流量作為噴管流動的關(guān)鍵參數(shù)之一，可在一定程度上衡量模型與實際的近似程度。為此將流量值與相關(guān)文獻(xiàn)實驗數(shù)據(jù)進(jìn)行對比。網(wǎng)格步長設(shè)定為-0.025 K，噴管喉部直徑為1.014 mm，其余初始條件設(shè)置如表4所示。

表4 流量驗證計算初始條件

流量計算值與實驗值對比如圖4所示。

圖4 噴管流動模型流量驗證Fig.4 Mass flow verification of the nozzle flow model

圖4中橫坐標(biāo)為噴管滯止?fàn)顟B(tài)下對比壓強(qiáng)：

以實驗數(shù)據(jù)作為基準(zhǔn)，流量計算值誤差如表5所示。

表5 流量計算誤差

根據(jù)圖4及表5可知：計算絕對誤差隨對比壓強(qiáng)增加而增加；相對誤差呈波動趨勢，始終保持在10%～20%的范圍內(nèi)。在實驗本身也可能存在測量儀器等多個誤差來源的前提下認(rèn)為計算誤差在可接受范圍內(nèi)。沿流動方向的比熵隨馬赫數(shù)變化如圖5所示。

圖5 沿流動方向比熵變化Fig.5 Specific entropy variation along the flow direction

由于本文關(guān)注的重點是噴管的臨界流量，比熵變化只考慮了由滯止?fàn)顟B(tài)加速至聲速的流動。根據(jù)圖5可知，在不同滯止?fàn)顟B(tài)下隨著二氧化碳加速流動，沿流向熵值呈現(xiàn)一定波動，但相對于滯止?fàn)顟B(tài)的最大變化量保持在5%以內(nèi)。因此認(rèn)為等熵假設(shè)可以滿足工程計算需求，噴管流動模型可行性和合理性得到驗證。

1.4 彈射內(nèi)彈道動力學(xué)模型

1.4.1 內(nèi)彈道模型

以下標(biāo)1表示高壓室，基于質(zhì)量守恒原理、能量守恒原理、真實氣體狀態(tài)方程和比焓定義式建立的高壓室動力學(xué)模型為

(28)

式中：為彈射初始時刻高壓室中工質(zhì)質(zhì)量；為高壓室壓強(qiáng)；為高壓室工質(zhì)定壓比熱容；為高壓室工質(zhì)定容比熱容；為高壓室體積；為高壓室工質(zhì)比體積；為高壓室溫度；為高壓室工質(zhì)比焓。

以下標(biāo)m表示包含彈托的導(dǎo)彈，下標(biāo)2表示低壓室，基于導(dǎo)彈運動學(xué)和動力學(xué)關(guān)系、質(zhì)量守恒原理、能量守恒原理、真實氣體狀態(tài)方程和比焓定義式建立的低壓室動力學(xué)模型為

(29)

式中：、分別為導(dǎo)彈位移和速度；為導(dǎo)彈質(zhì)量；為低壓室有效推力面積；為低壓室壓強(qiáng)；為大氣壓；為重力加速度；為低壓室體積；為低壓室工質(zhì)質(zhì)量；為低壓室溫度；為低壓室工質(zhì)比體積；為低壓室工質(zhì)定容比熱容；為低壓室工質(zhì)比焓。

142 內(nèi)彈道模型驗證

將基于超臨界空氣的垂直發(fā)射裝置內(nèi)彈道計算結(jié)果與文獻(xiàn)數(shù)據(jù)進(jìn)行了對比，干空氣熱物性參數(shù)如表6所示。

表6 干空氣熱物性參數(shù)

用于內(nèi)彈道模型驗證的仿真參數(shù)如表7所示。

表7 內(nèi)彈道模型驗證仿真參數(shù)

內(nèi)彈道模型驗證結(jié)果如圖6所示。

圖6 基于空氣的內(nèi)彈道計算結(jié)果與文獻(xiàn)[2]數(shù)據(jù)對比Fig.6 Comparison of the calculated interior ballistics results based on air with the data from Ref.[2]

根據(jù)圖6，基于前文內(nèi)彈道模型計算的導(dǎo)彈加速度、速度、位移等參數(shù)與文獻(xiàn)[2]數(shù)據(jù)相比，吻合較好。計算值與文獻(xiàn)[2]數(shù)值的變化趨勢一致，具體數(shù)值存在一定差別。以文獻(xiàn)[2]數(shù)值作為基準(zhǔn)，則計算的加速度最大絕對差值為5.178 m/s，對應(yīng)的相對差值為7.56%；速度最大絕對差值1.622 m/s，對應(yīng)的相對差值為13.74%；位移最大絕對差值0.672 m，對應(yīng)的相對差值為5.47%。

綜合上述分析可知，考慮到文獻(xiàn)[2]中使用的建模方法與本文有所區(qū)別，認(rèn)為計算結(jié)果差值在可接受范圍之內(nèi)，本文使用的內(nèi)彈道模型合理性和可行性得到驗證。

2 內(nèi)彈道動力學(xué)仿真與優(yōu)化

2.1 內(nèi)彈道仿真計算

2.1.1 時間步長對計算影響

為研究時間步長d對內(nèi)彈道模型仿真結(jié)果的影響，工質(zhì)選定為超臨界二氧化碳，在2×10～1×10s的范圍內(nèi)，以2×10s為間隔設(shè)計5個計算工況。設(shè)置噴管計算步長為-0.037 K，噴管喉部直徑在0～0.4 s為85 mm，在0.4～0.6 s為120.2 mm，0.6 s之后為147.2 mm。計算初始條件及導(dǎo)彈參數(shù)如表8所示，取9.81 m/s。

表8 時間步長對比計算參數(shù)

基于超臨界二氧化碳的內(nèi)彈道計算結(jié)果如圖7所示。

圖7 不同時間步長計算結(jié)果對比Fig.7 Calculation results using different time steps

根據(jù)圖7(a)可知，彈射開始后高壓室壓強(qiáng)由于工質(zhì)流出而持續(xù)降低，低壓室壓強(qiáng)先升高后降低，峰值壓強(qiáng)約為0.568 MPa。噴管喉部直徑的變化有效減緩了低壓室壓強(qiáng)的下降。結(jié)合加速度變化曲線可知，改變噴管喉部橫截面積可實現(xiàn)對彈射過程的控制。根據(jù)圖7(b)，高壓室溫度由初始373 K降至約295 K，低壓室溫度先增大至約340 K后持續(xù)降至約225 K。由(29)式可知，在0.08 s之前高壓室與低壓室比焓差值相對較小，低壓室溫度對時間的導(dǎo)數(shù)為正值，使得低壓室溫度升高。根據(jù)圖7(c)，導(dǎo)彈加速度峰值約5.05，未超過過載限制。根據(jù)圖7(d)和圖7(e)，導(dǎo)彈速度和位移隨時間變化持續(xù)增加直至彈射過程終止。終止時刻的導(dǎo)彈運動學(xué)參數(shù)如表9所示。

表9 彈射終止時刻導(dǎo)彈運動學(xué)參數(shù)

結(jié)合圖7及表9可知，不同算例對應(yīng)的高、低壓室狀態(tài)量及導(dǎo)彈運動學(xué)參數(shù)基本一致，表明內(nèi)彈道模型對計算步長的變化不敏感，后續(xù)可在綜合考慮求解精度及時間的基礎(chǔ)上設(shè)置計算步長。

2.1.2 真實氣體效應(yīng)對計算影響

為研究真實氣體效應(yīng)對模型影響，對比了分別采用真實氣體和理想氣體模型的內(nèi)彈道計算結(jié)果。仿真參數(shù)及內(nèi)彈道方程組初始條件與2.1.1節(jié)中設(shè)置基本一致，計算步長取5×10s。計算結(jié)果如圖8所示。

圖8 理想氣體和真實氣體模型計算結(jié)果對比Fig.8 Comparison of calculation results for ideal gas and real gas models

由圖8(a)可知，高壓室壓強(qiáng)持續(xù)降低；低壓室壓強(qiáng)在彈射初期上升，之后由于工質(zhì)流入的增壓效應(yīng)無法抵消體積增大的影響而開始降低?；诶硐霘怏w模型計算的低壓室壓強(qiáng)峰值約0.626 MPa，明顯大于真實氣體模型計算結(jié)果0.568 MPa，相對差值約10.35%。根據(jù)(29)式即低壓室內(nèi)彈道方程組分析原因，高、低壓室比焓差值(-)一項對低壓室壓強(qiáng)、溫度和比焓變化率均具有重要影響。根據(jù)理想氣體假設(shè)，工質(zhì)比焓隨溫度增加而增加。在彈射初期，高壓室溫度高于低壓室，導(dǎo)致由理想模型計算的高壓室比焓顯著大于低壓室比焓。較大的比焓差值使得低壓室壓強(qiáng)峰值有較大幅度增加。

由圖8(b)可知，彈射初始階段高壓室、低壓室間較大的比焓差值也使得基于理想氣體模型計算的低壓室溫度峰值與真實氣體模型計算結(jié)果相比，增加約100 K左右。同時，綜合圖8(a)和圖8(b)可知，真實氣體效應(yīng)對高壓室壓強(qiáng)、溫度影響不明顯。

由圖8(c)可知，真實氣體模型流量計算值大于理想氣體計算結(jié)果，表明真實氣體效應(yīng)對噴管內(nèi)工質(zhì)流動存在影響。分析原因，主要是在相同的滯止壓強(qiáng)和溫度下，基于理想氣體模型計算的工質(zhì)密度顯著低于真實氣體模型計算值，同時不同模型計算的無量綱流動系數(shù)也存在差別。

由圖8(d)可知，根據(jù)理想氣體假設(shè)計算的導(dǎo)彈加速度在0～0.4 s之間明顯偏大，峰值加速度為56.978 m/s。根據(jù)真實氣體模型計算的峰值加速度為49.505 m/s，區(qū)別主要是由于理想氣體模型計算的較高低壓室壓強(qiáng)造成的。加速度對比也表明理想氣體模型下導(dǎo)彈運動存在超過過載限制的可能。彈射后期理想氣體模型加速度下降速率與真實氣體模型相比，有較小幅度增加。

綜上所述，工質(zhì)的真實氣體效應(yīng)對內(nèi)彈道模型特別是低壓室部分影響較大，建模和計算需要充分考慮相關(guān)因素。

2.2 內(nèi)彈道優(yōu)化設(shè)計

2.2.1 內(nèi)彈道優(yōu)化模型

為提升彈射整體效能，在仿真基礎(chǔ)上設(shè)計了內(nèi)彈道優(yōu)化模型。由于單一評價指標(biāo)無法充分體現(xiàn)不同優(yōu)化方案的差別，引入多目標(biāo)優(yōu)化思想。

2.2.1.1 優(yōu)化目標(biāo)設(shè)定

優(yōu)化目標(biāo)一般要結(jié)合問題的具體背景確定。對于彈射過程，首要關(guān)注點在于能否實現(xiàn)給定的彈射指標(biāo)即出筒速度。但由于并未嚴(yán)格限定出筒速度和設(shè)計值之間的偏差，出筒速度可以作為約束條件而非優(yōu)化目標(biāo)引入優(yōu)化模型中。在此基礎(chǔ)上彈射裝置的整體效能，亦即為了實現(xiàn)預(yù)期目標(biāo)所使用的工質(zhì)和工質(zhì)能量轉(zhuǎn)化率成為關(guān)注重點。根據(jù)這一思想設(shè)計如下優(yōu)化目標(biāo)：

1)基于工質(zhì)質(zhì)量的目標(biāo)函數(shù)。彈射裝置使用的超臨界二氧化碳質(zhì)量一方面影響整體使用成本，另一方面決定了彈射裝置體積。較大的體積將顯著增加武器系統(tǒng)成本和受攻擊風(fēng)險。所需工質(zhì)質(zhì)量應(yīng)當(dāng)在滿足彈射指標(biāo)的情況下盡可能較少，同時為了使各目標(biāo)函數(shù)的數(shù)量級保持一致，與工質(zhì)質(zhì)量相關(guān)的優(yōu)化目標(biāo)可表示為

(30)

2)基于能量利用率的目標(biāo)函數(shù)。工質(zhì)的能量利用效率可用于評估彈射裝置效能。給定工質(zhì)能量利用率為

式中：為導(dǎo)彈離筒時刻機(jī)械能；u為初始時刻高壓室內(nèi)工質(zhì)內(nèi)能?？偟哪芰坷寐蕬?yīng)當(dāng)保持在相對較高的水平：

=min(1-)

(31)

2212 設(shè)計變量選取

設(shè)計變量的選取參考相關(guān)文獻(xiàn)工作。高壓室和低壓室初始狀態(tài)直接決定了彈射動力學(xué)特性，影響彈射效能。例如高壓室初始壓強(qiáng)、溫度和體積共同決定了使用工質(zhì)質(zhì)量；低壓室初始體積影響低壓室內(nèi)壓強(qiáng)變化。因此選定高壓室初始壓強(qiáng)、高壓室初始溫度、高壓室體積、低壓室初始體積等4個變量作為設(shè)計變量。約束條件如下：

1)設(shè)計變量上下限約束。設(shè)計變量的取值范圍應(yīng)當(dāng)結(jié)合問題物理背景進(jìn)行限制：

(32)

式中：、上下限主要根據(jù)S-R-K方程精度設(shè)置；、上下限根據(jù)文獻(xiàn)[1]的結(jié)論和試算確定。

2)噴管臨界流動約束。根據(jù)內(nèi)彈道模型假設(shè)，噴管流動始終處于臨界狀態(tài)。要求噴管出口和入口、即低壓室和高壓室壓強(qiáng)之比不得大于臨界值：

(33)

式中：為臨界對比壓比，等于噴管喉部壓強(qiáng)與噴管入口處滯止壓強(qiáng)之比。實際計算時該臨界值由基于真實氣體效應(yīng)的一維等熵噴管流動模型確定。

3)低壓室工質(zhì)三相點約束。在彈射過程中，低壓室工質(zhì)溫度存在下降階段，可能降至工質(zhì)三相點溫度以下。此時若工質(zhì)壓強(qiáng)大于三相點壓強(qiáng)，則相態(tài)為固態(tài)，可能造成彈射失敗。該約束可表示為

<,<

(34)

式中：、分別為工質(zhì)三相點的壓強(qiáng)、溫度。

4)導(dǎo)彈最大過載約束。彈射過程中由于彈體結(jié)構(gòu)和彈上設(shè)備的要求，導(dǎo)彈加速度不得超過許用值：

≤

(35)

5)導(dǎo)彈出筒速度約束。導(dǎo)彈出筒速度應(yīng)大于等于設(shè)計值：

≥

(36)

(30)式～(36)式即為內(nèi)彈道優(yōu)化模型。

2.2.2 基于多目標(biāo)遺傳算法的優(yōu)化模型求解

基于多目標(biāo)遺傳算法——帶精英策略的非支配排序的遺傳算法對模型進(jìn)行求解。仿真參數(shù)和彈射設(shè)計指標(biāo)主要根據(jù)表8參數(shù)設(shè)置，內(nèi)彈道模型計算步長為1×10s，低壓室初始時刻溫度為283 K。為了提高導(dǎo)彈出筒時刻機(jī)動性，出筒速度設(shè)計值調(diào)整為30 m/s。遺傳算法相關(guān)參數(shù)設(shè)置為：種群數(shù)量50，迭代次數(shù)最大取120，交叉概率0.8。

隨著迭代次數(shù)增加，Pareto最優(yōu)解集對應(yīng)的目標(biāo)函數(shù)值變化如圖9所示。

圖9 Pareto最優(yōu)解集對應(yīng)目標(biāo)函數(shù)變化Fig.9 Changes in the objective function corresponding to the Pareto optimal set

根據(jù)圖9(a)可知，單個目標(biāo)函數(shù)的最小值均隨迭代次數(shù)增加而呈現(xiàn)下降趨勢。相對迭代次數(shù)為60時，迭代次數(shù)120的情況下目標(biāo)函數(shù)下降0.000 1%，目標(biāo)函數(shù)下降0.001 9%，表明優(yōu)化目標(biāo)對迭代次數(shù)的變化更為敏感。根據(jù)圖9(b)，隨迭代次數(shù)增加，目標(biāo)函數(shù)Pareto前沿在目標(biāo)空間內(nèi)沿函數(shù)值減小方向的推進(jìn)比較明顯，表明優(yōu)化取得一定效果。同時由于優(yōu)化目標(biāo)的改進(jìn)量相對較小，認(rèn)為在現(xiàn)有條件下迭代次數(shù)設(shè)置在120的范圍內(nèi)是較為合理的。

120次迭代的Pareto解及對應(yīng)目標(biāo)函數(shù)值與特定基準(zhǔn)工況對比如表10所示。

表10 迭代120次條件下優(yōu)化結(jié)果對比

根據(jù)表10可知：與優(yōu)化前相比，優(yōu)化后的高壓室壓強(qiáng)增大21.09%，溫度近似相等；高壓室體積減小46.33%，低壓室初始體積增加約44.67%。由于并未嚴(yán)格限定彈射裝置整體的結(jié)構(gòu)尺寸和安全性指標(biāo)，其在當(dāng)前問題中均屬于次要因素，因此認(rèn)為優(yōu)化解相對合理。

目標(biāo)函數(shù)方面，基礎(chǔ)解所需的工質(zhì)質(zhì)量約為1 226 kg，優(yōu)化后所需的工質(zhì)質(zhì)量減小至約831 kg，與基礎(chǔ)解相比減少32.22%。能量利用率由2.10%增加至3.14%，增加49.52%。優(yōu)化前后導(dǎo)彈運動參數(shù)曲線如圖10所示。

圖10 優(yōu)化前后導(dǎo)彈運動學(xué)參數(shù)對比Fig.10 Comparison of missile kinematic parameters before and after optimization

根據(jù)圖10可知，基準(zhǔn)解所需彈射時間約0.965 s，優(yōu)化后彈射時間增加至約1.011 s，相對優(yōu)化前增加約4.76%。綜合實際彈射過程認(rèn)為，彈射時間增加所造成的影響可以近似忽略。優(yōu)化后導(dǎo)彈加速度曲線在上升階段更為平緩，加速度峰值到達(dá)時間延后。優(yōu)化后加速度峰值約43.892 m/s，與優(yōu)化前峰值43.443 m/s相比近似相等。優(yōu)化后加速度曲線的相對平緩使得出筒時間有所延長，但如前所述對系統(tǒng)反應(yīng)能力的影響可忽略不計。

3 結(jié)論

本文針對使用超臨界二氧化碳作為工作介質(zhì)的導(dǎo)彈彈射裝置，基于S-R-K狀態(tài)方程構(gòu)建了內(nèi)彈道動力學(xué)模型。在研究時間步長和真實氣體效應(yīng)對模型仿真結(jié)果影響的基礎(chǔ)上，設(shè)計了多目標(biāo)優(yōu)化模型并進(jìn)行求解。得到如下主要結(jié)論：

1)內(nèi)彈道模型精度通過實驗數(shù)據(jù)和相關(guān)文獻(xiàn)仿真結(jié)果得到驗證，表明本文建立的內(nèi)彈道模型較為合理，模型采用的數(shù)值求解方法是有效的。

2)時間步長對計算結(jié)果影響不明顯，表明可在綜合精度和時間成本等因素的基礎(chǔ)上選擇合適的時間步長進(jìn)行求解。

3)對比了采用理想氣體模型和真實氣體模型的計算結(jié)果，理想氣體假設(shè)下低壓室溫度和壓強(qiáng)在彈射過程中明顯偏高，表明真實氣體效應(yīng)對計算存在比較顯著的影響。

4)基于遺傳算法的優(yōu)化結(jié)果與參考值相比明顯改進(jìn)，工質(zhì)用量減少32.22%，能量利用率提升49.52%。表明遺傳算法對基于超臨界二氧化碳的彈射器內(nèi)彈道多目標(biāo)優(yōu)化問題具有一定適用性，優(yōu)化后的彈射系統(tǒng)多項指標(biāo)均有改進(jìn)。

后續(xù)可在引入多相共存及非平衡效應(yīng)的基礎(chǔ)上建立更為精細(xì)的內(nèi)彈道模型，使計算結(jié)果更加準(zhǔn)確。

[1] 譚大成.彈射內(nèi)彈道學(xué)[M].北京:北京理工大學(xué)出版社，2015.

TAN D C.Interior ballistics of catapult[M].Beijing:Beijing Institute of Technology Press，2015.(in Chinese)

[2] REN J，YANG F B，MA D W，et al.Pneumatic performance study of a high pressure ejection device based on real specific energy and specific enthalpy[J].Entropy，2014，16(9):4801-4817.

[3] 姚琳，馬大為，馬吳寧，等.兩級提拉式單側(cè)彈射裝置內(nèi)彈道建模與優(yōu)化[J].兵工學(xué)報，2017，38(3):466-475.

YAO L，MA D W，MA W N，et al.Interior ballistics modeling and optimization of one-side ejection device with two-step cylinder[J].Acta Armamentarii，2017，38(3):466-475. (in Chinese)

[4] YAO H X，WEI X Y，YE H. Supercritical carbon dioxide as a new working medium for pneumatic launch: a theoretical study[J]. Defence Technology，2020，17(4):1296-1306.

[5] LETTIERI C，PAXSON D，SPAKOVSZKY Z，et al. Characterization of nonequilibrium condensation of supercritical carbon dioxide in a de laval nozzle[J].Journal of Engineering for Gas Turbines and Power，2018，140(4): 041701.

[6] SOAVE G.Equilibrium constants from a modified Redlich-Kwong equation of state[J].Chemical Engineering Science，1972，27(6): 1197-1203.

[7] SIRIGNANO W A.Compressible flow at high pressure with linear equation of state[J].Journal of Fluid Mechanics，2018，843:244-292.

[8] GRABOSKI M S，DAUBERT T E.A modified Soave equation of state for phase equilibrium calculations.3. Systems containing hydrogen[J].Industrial & Engineering Chemistry Process Design and Development，1979，18(2):300-306.

[9] BOBER W，CHOW W L.Nonideal isentropic gas flow through converging-diverging nozzles[J].Journal of Fluids Engineering，1990，112(4):455-460.

[10] BLUNDELL S J，BLUNDELL K M.Concepts in thermal physics[M].Oxford，UK: Oxford University Press，2009.

[11] PASSUT C A，DANNER R P.Correlation of ideal gas enthalpy，heat capacity and entropy[J].Industrial & Engineering Chemistry Process Design and Development，1972，11(4): 543-546.

[12] NIST.NIST chemistry WebBook，SRD69[DB/OL].Gaithersburg:National Institute of Standards and Technology，2018[2021-05-10].https://webbook.nist.gov/chemistry/fluid/.

[13] BEUNE A.Analysis of high-pressure safety valves[D]. Eindhoven，the Nethlands:Eindhoven University of Technology，2009.

[14] MORAN M J，SHAPIRO H N.Fundamentals of engineering thermodynamics[M].5th edit.Hoboken，NJ, US:John Wiley & Sons，2006.

[15] EDLEBECK J，NELLIS G F，KLEIN S A，et al.Measurements of the flow of supercritical carbon dioxide through short orifices[J].The Journal of Supercritical Fluids，2014，88:17-25.

[16] 馬營，楊臻，趙蔚楠，等.內(nèi)彈道多目標(biāo)優(yōu)化算法對比分析[J].兵工自動化，2021，40(3): 81-85.

MA Y，YANG Z，ZHAO W N，et al.Comparative analysis of multiple objective optimization algorithms for interior ballistics[J].Ordnance Industry Automation，2021，40(3):81-85. (in Chinese)