基于超臨界二氧化碳的導(dǎo)彈彈射裝置內(nèi)彈道建模與優(yōu)化

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

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

0 引言

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

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

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

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

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

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

1.1 模型假設(shè)

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

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

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

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

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

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

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

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

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

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

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

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

理論建模中一般通過(guò)狀態(tài)方程引入真實(shí)氣體效應(yīng)。S-R-K狀態(tài)方程是Soave基于Redlich-Kwong方程提出的一種立方型狀態(tài)方程，以在多個(gè)應(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)量計(jì)算表達(dá)式。計(jì)算公式為

(3)

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

(4)

(5)

(6)

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

(7)

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

(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)式，即可計(jì)算特定壓強(qiáng)和溫度下工質(zhì)定壓比熱容。

比焓可表示為

(11)

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

(12)

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

(13)

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

(14)

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

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

(15)

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

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

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

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

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

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

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

圖2 基于S-R-K方程的二氧化碳熱物性參數(shù)計(jì)算值與NIST數(shù)據(jù)對(duì)比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)時(shí)，S-R-K狀態(tài)方程預(yù)測(cè)的二氧化碳物性參數(shù)計(jì)算誤差相對(duì)較小，可以滿足后續(xù)內(nèi)彈道計(jì)算需求。

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

1.3.1 一維等熵噴管流動(dòng)模型

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

(16)

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

(17)

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

(18)

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

132 一維等熵噴管流動(dòng)的迭代解法

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

(19)

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

(20)

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

(21)

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

(22)

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

=-1+d

(23)

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

(24)

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

(25)

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

(26)

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

(27)

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

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

133 噴管流動(dòng)模型驗(yàn)證

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

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

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

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

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

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

表4 流量驗(yàn)證計(jì)算初始條件

流量計(jì)算值與實(shí)驗(yàn)值對(duì)比如圖4所示。

圖4 噴管流動(dòng)模型流量驗(yàn)證Fig.4 Mass flow verification of the nozzle flow model

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

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

表5 流量計(jì)算誤差

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

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

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

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

1.4.1 內(nèi)彈道模型

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

(28)

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

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

(29)

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

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

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

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

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

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

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

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

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

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

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

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

2.1.1 時(shí)間步長(zhǎng)對(duì)計(jì)算影響

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

表8 時(shí)間步長(zhǎng)對(duì)比計(jì)算參數(shù)

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

圖7 不同時(shí)間步長(zhǎng)計(jì)算結(jié)果對(duì)比Fig.7 Calculation results using different time steps

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

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

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

2.1.2 真實(shí)氣體效應(yīng)對(duì)計(jì)算影響

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

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

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

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

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

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

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

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

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

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

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

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

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

(30)

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

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

=min(1-)

(31)

2212 設(shè)計(jì)變量選取

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

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

(32)

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

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

(33)

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

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

<,<

(34)

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

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

≤

(35)

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

≥

(36)

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

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

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

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

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

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

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

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

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

目標(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)彈運(yùn)動(dòng)參數(shù)曲線如圖10所示。

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

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

3 結(jié)論

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

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

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

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

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

后續(xù)可在引入多相共存及非平衡效應(yīng)的基礎(chǔ)上建立更為精細(xì)的內(nèi)彈道模型，使計(jì)算結(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] 姚琳，馬大為，馬吳寧，等.兩級(jí)提拉式單側(cè)彈射裝置內(nèi)彈道建模與優(yōu)化[J].兵工學(xué)報(bào)，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] 馬營(yíng)，楊臻，趙蔚楠，等.內(nèi)彈道多目標(biāo)優(yōu)化算法對(duì)比分析[J].兵工自動(dòng)化，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)