孫曉暉，陳志華，韓珺禮，2，薛大文

（1.南京理工大學(xué) 瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室，南京210094；2.北京機(jī)電研究所，北京100012）

沖壓加速器是由美國(guó) HERTZBERY[1]提出的一種高超聲速發(fā)射技術(shù)，它利用化學(xué)能加速?gòu)椡柚粮叱曀?，其推力穩(wěn)定，過載低，還克服了傳統(tǒng)火炮彈丸加速時(shí)彈底壓力減小的缺點(diǎn)，具有高效的化學(xué)能與動(dòng)能轉(zhuǎn)化率，費(fèi)用低廉，因而在超高聲速氣動(dòng)物理和氣動(dòng)熱力學(xué)研究、航天和軍事飛行器發(fā)射及超燃沖壓發(fā)動(dòng)機(jī)機(jī)理研究與空間輸運(yùn)等方面具有很好的應(yīng)用前景.

BRUCKNER等[2]實(shí)驗(yàn)研究發(fā)現(xiàn)，提高填充壓力并采用含能更高的推進(jìn)劑將使彈丸獲得更大的速度增量；HERTZBERG等[3]研究了推力隨彈丸速度的變化情況，發(fā)現(xiàn)在熱壅塞模式下實(shí)驗(yàn)結(jié)果與理論模型相符，當(dāng)彈丸速度接近推進(jìn)劑CJ爆轟速度時(shí)，彈丸加速情況則優(yōu)于理論預(yù)期；BENGHERBIA等[4，5]通過五步化學(xué)反應(yīng)模式，數(shù)值模擬了沖壓加速器亞爆轟模態(tài)，預(yù)測(cè)了推力隨馬赫數(shù)的變化趨勢(shì)，并與實(shí)驗(yàn)情況基本相符.

國(guó)內(nèi)學(xué)者同樣對(duì)沖壓加速器進(jìn)行了大量研究，如柳森等[6，7]利用氣動(dòng)中心研制的37mm口徑?jīng)_壓加速裝置將彈丸加速至1 760m/s以上，提供了熱壅塞模式管壁的壓力分布與彈丸速度變化情況等實(shí)驗(yàn)數(shù)據(jù)，為沖壓加速器的數(shù)值研究提供了參考，并開展了冷態(tài)實(shí)驗(yàn)加速管內(nèi)彈丸繞流的數(shù)值模擬.陳堅(jiān)強(qiáng)、張涵信等[8]利用ENN格式，數(shù)值摸擬了 H2/Air斜爆轟波沖壓加速器繞流，指出粘性效應(yīng)對(duì)燃燒過程具有重要的影響，燃燒首先發(fā)生在邊界層中，并依賴于來流條件，可發(fā)生穩(wěn)定和不穩(wěn)定2種燃燒過程.張國(guó)強(qiáng)、翁春生等[9，10]通過數(shù)值研究表明預(yù)混可燃?xì)怏w在彈后點(diǎn)燃，并在彈后某一位置形成全管面燃燒.

準(zhǔn)確、精細(xì)的流場(chǎng)結(jié)構(gòu)對(duì)系統(tǒng)設(shè)計(jì)、彈丸形狀與加速機(jī)理等可提供重要參考，然而受計(jì)算格式與資源的影響，相關(guān)沖壓加速器冷態(tài)實(shí)驗(yàn)數(shù)值計(jì)算結(jié)果無法有效地揭示流場(chǎng)中復(fù)雜激波流場(chǎng)細(xì)節(jié)及重要參數(shù)分布與變化情況.基于此，本文利用高精度加權(quán)本質(zhì)無震蕩（Weight Essentially Non－Oscillatory，WENO）計(jì)算格式對(duì)二維可壓歐拉方程進(jìn)行求解，結(jié)合自適應(yīng)網(wǎng)格加密技術(shù)（Adaptive Mesh Refinement，AMR）與沉浸邊界法（Immersed Boundary Method，IBM）對(duì)方形沖壓加速器[11]冷態(tài)實(shí)驗(yàn)進(jìn)行數(shù)值模擬，并與相應(yīng)的實(shí)驗(yàn)流場(chǎng)紋影圖進(jìn)行對(duì)比，驗(yàn)證了數(shù)值方法對(duì)激波流場(chǎng)細(xì)節(jié)捕捉的可行性.同時(shí)計(jì)算結(jié)果清晰地展示了流場(chǎng)的復(fù)雜波系結(jié)構(gòu)與重要物理參數(shù)的變化情況，對(duì)沖壓加速器冷態(tài)實(shí)驗(yàn)相關(guān)研究具有一定參考價(jià)值.

1 計(jì)算方法

本文采用二維可壓歐拉方程：

式中，

式中，ρ為氣體密度；u、v分別為x與y方向的速度分量；p為壓力；T為溫度；E為單位質(zhì)量氣體的總能量，E＝[p/ρ（r－1）]＋0.5（u2＋v2），r為理想氣體絕熱指數(shù).理想氣體的狀態(tài)方程為p＝ρRT，R為理想氣體常數(shù).

WENO 格式[12]是在 ENO（Essentially Non－oscillatory）格式的基礎(chǔ)上，通過引入變化的加權(quán)因子對(duì)不同模板的光滑程度進(jìn)行優(yōu)化組合，以提高光滑區(qū)域的計(jì)算精度，并可在間斷區(qū)域保持ENO格式的分辨率.因此，WENO格式具有很好的有效性、通量的光滑性與收斂解的穩(wěn)定性，并保持對(duì)激波類的移動(dòng)間斷和接觸界面的捕捉能力.

本文采用笛卡爾網(wǎng)格，它具有網(wǎng)格建立簡(jiǎn)單、快速、數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單、網(wǎng)格自適應(yīng)容易的特點(diǎn)，可縮短流場(chǎng)計(jì)算時(shí)間，且網(wǎng)格生成自動(dòng)化簡(jiǎn)單.AMR技術(shù)可根據(jù)計(jì)算物理量梯度變化自動(dòng)調(diào)整當(dāng)?shù)鼐W(wǎng)格大小，以保證對(duì)流場(chǎng)中參數(shù)變化大的區(qū)域的計(jì)算精度，從而盡可能用較小的總網(wǎng)格數(shù)對(duì)所求解計(jì)算域的物理問題給出較高精度的數(shù)值解，并可在同等計(jì)算精度的情況下大大減少計(jì)算量[13].笛卡爾網(wǎng)格只能用在幾何結(jié)構(gòu)簡(jiǎn)單的計(jì)算域中，為解決此限制，結(jié)合IBM方法對(duì)外形較為復(fù)雜的邊界進(jìn)行處理.IBM[14]的基本思想是對(duì)研究流場(chǎng)進(jìn)行拓展，使流場(chǎng)邊界完全包含在拓展流場(chǎng)中，原流場(chǎng)邊界變?yōu)樾铝鲌?chǎng)內(nèi)點(diǎn)，因此無需考慮邊界形狀，可對(duì)復(fù)雜邊界進(jìn)行數(shù)值計(jì)算.但I(xiàn)BM需保證拓展流場(chǎng)的合理性，特別是對(duì)原邊界條件的滿足.相關(guān)技術(shù)的詳細(xì)描述參見文獻(xiàn)[14].

為了驗(yàn)證以上數(shù)值方法的可行性，以激波繞過三角楔過程的激波流場(chǎng)（SCHARDIN問題）為算例進(jìn)行數(shù)值驗(yàn)證.實(shí)驗(yàn)陰影[15]與本文的數(shù)值結(jié)果比較如圖1所示.可知兩者幾乎完全相符.激波繞射過程中，整個(gè)流場(chǎng)基本呈軸對(duì)稱，入射激波繞過三角楔并與其相互作用，在三角楔尾部的上、下角處分別形成順時(shí)針和逆時(shí)針的大渦，同時(shí)上、下反射激波分別與大渦相互作用.數(shù)值結(jié)果清晰地展示了流場(chǎng)中大渦和渦串結(jié)構(gòu)，及其與激波的相互作用，表明了以上方法對(duì)激波流場(chǎng)的高分辨率.

圖1 SCHARDIN問題的實(shí)驗(yàn)與本文計(jì)算結(jié)果比較

2 結(jié)果與討論

本文彈丸尺寸及其附近網(wǎng)格分布如圖2所示.可知，彈丸與管壁構(gòu)建成一個(gè)類似噴管的結(jié)構(gòu)，其最小截面處稱為管道喉部.彈丸大小參考文獻(xiàn)[11]，具體標(biāo)注如圖2（a）所示，其彈丸頭部固定在（0，0）處.

取計(jì)算域大小為0.17m×0.023m，其x方向與y方向坐標(biāo)范圍分別為（－0.020，0.15）與（－0.011 5，0.011 5）.設(shè)管壁與彈丸表面為絕熱反射壁面，右端為開口條件.左端初始來流為空氣，具體參數(shù)：壓力p＝0.4 MPa，溫度T＝303K，來流馬赫數(shù)Ma∞＝3.5.初始笛卡爾網(wǎng)格大小為Δx＝0.125mm，Δy＝0.125 mm.網(wǎng)格自適應(yīng)加密層數(shù)為2，加密因子為2.

圖2 初始時(shí)刻彈丸周圍的網(wǎng)格分布與局部放大圖

圖3為文獻(xiàn)[11]的冷態(tài)實(shí)驗(yàn)紋影與流場(chǎng)結(jié)構(gòu)及本文數(shù)值計(jì)算紋影的比較.實(shí)驗(yàn)條件：管中為CH4＋2O2＋4.8CO2混合氣體，其溫度為303 K，壓力為0.4 MPa，混合氣體聲速為300m/s，彈丸飛行速度為1 210m/s，因而彈丸飛行馬赫數(shù)為Ma∞＝4.04.可知，本文數(shù)值計(jì)算紋影與實(shí)驗(yàn)紋影基本相符.圖中斜激波、膨脹波系以及激波在壁面的反射清晰可見.相比實(shí)驗(yàn)圖片，數(shù)值計(jì)算的波系結(jié)構(gòu)更加清晰，并能展示彈丸下游流場(chǎng)中的旋渦與激波的相互作用以及激波相交所形成的復(fù)雜波系結(jié)構(gòu).其中，AB稱為前楔激波、BD為前楔反射激波、DG為后楔反射激波、IJ為彈底激波.

圖3 實(shí)驗(yàn)的流場(chǎng)紋影圖和流場(chǎng)結(jié)構(gòu)示意圖及本文數(shù)值計(jì)算的流場(chǎng)紋影圖

圖4為不同時(shí)間段，管內(nèi)流場(chǎng)發(fā)展變化的密度等值線分布.超聲速氣流在彈丸頭部形成對(duì)稱斜激波，同時(shí)受頭部楔面（前楔）氣流初始運(yùn)動(dòng)的影響，表面有法向正激波，并使斜激波陣面稍微向外彎曲，此時(shí)激波與管壁碰撞反射角度較小，彈丸肩部形成典型普朗特－邁耶膨脹.彈丸尾部楔面（后楔）及尾部初始流場(chǎng)結(jié)構(gòu)剛形成，彈丸附近流場(chǎng)的激波強(qiáng)度都較弱，如圖4（a）所示.隨著運(yùn)動(dòng)發(fā)展，前楔斜激波強(qiáng)度得到提高，傾角開始固定，同時(shí)波陣面駐定并變?yōu)橹本€，且與管壁的碰撞反射點(diǎn)前移，此時(shí)前楔斜激波與管壁的反射角變大且穩(wěn)定在20°左右（圖4（b）～圖4（e）），反射使波后氣體壓力、溫度急劇升高，但在彈丸肩部的膨脹波作用下，波后氣體的溫度和壓力迅速恢復(fù)到正常水平，同時(shí)在彈肩后方形成局部低壓區(qū).

圖4 管內(nèi)流場(chǎng)發(fā)展過程的密度等值線圖（單位：kg/m3）

前楔反射激波剛形成時(shí)強(qiáng)度較弱，波陣面呈彎曲狀，隨氣流不斷向后楔表面運(yùn)動(dòng)，在t＝0.04ms時(shí)，達(dá)到后楔表面并發(fā)生反射，形成后楔反射激波.此時(shí)，前楔反射激波陣面幾乎成直線，在彈肩后方低壓區(qū)的作用下，波陣面底部發(fā)生彎曲，如圖4（b）所示.隨著流場(chǎng)的發(fā)展，最終在t＝0.06ms時(shí)，前楔反射激波傾角固定，反射點(diǎn)前移并穩(wěn)定在后楔表面.后楔反射激波繼續(xù)向下游傳播，同時(shí)向管道壁面偏移，在管道壁面發(fā)生反射，反射點(diǎn)最終固定在x＝0.075m處（圖4（e））.

氣流在繞經(jīng)彈丸尾部時(shí)同樣發(fā)生膨脹，在彈丸底部上下兩側(cè)形成低壓區(qū)，膨脹扇在管道軸線處發(fā)生碰撞并形成兩道較強(qiáng)的激波，稱為彈底激波.開始時(shí)，彈底激波的強(qiáng)度較弱，由于復(fù)雜波系的影響，波陣面十分粗糙.隨著氣體繼續(xù)向下游運(yùn)動(dòng)，彈底激波不斷強(qiáng)化，波陣面逐漸平滑，并在管道壁面形成反射，反射激波與反射點(diǎn)不斷向上游移動(dòng).在t＝0.24ms時(shí)，彈底激波基本穩(wěn)定，并在管道壁面發(fā)生兩次反射.

由于斜壓效應(yīng)，在彈丸底部上、下兩側(cè)形成兩道旋渦，旋渦上下相間隨機(jī)脫落，形成渦街，影響下游流場(chǎng).渦街在向下游傳播的過程中，不斷向兩側(cè)擴(kuò)散，最多至沖壓加速管徑的2/3，強(qiáng)度也不斷下降，但在彈底激波及其反射波的作用下得到強(qiáng)化與會(huì)聚，如圖4（d）～圖4（e）.各種波系與渦街相互作用，使得流場(chǎng)結(jié)構(gòu)極其復(fù)雜.

圖5為管內(nèi)流場(chǎng)溫度等勢(shì)分布圖，由于采用無粘滑移壁面，彈丸表面溫度不高，流場(chǎng)溫度極值（T≥1 150K）出現(xiàn)在彈丸底部，起到火焰穩(wěn)定器的作用，彈丸底部的渦街也保持較高的溫度，這有利于燃燒的擴(kuò)散與傳播.前楔斜激波在管道壁面與彈丸表面的反射使其波后的溫度升高，分別達(dá)到T≥498K、T≥465K，可以預(yù)熱甚至點(diǎn)燃預(yù)混可燃?xì)?在彈丸下游，后楔反射激波與彈底激波的兩次反射作用，使得管壁附近的溫度得到較大提升（T≥480K），可以促進(jìn)彈底燃燒擴(kuò)散至全管，對(duì)于沖壓加速器熱壅塞燃燒擴(kuò)散到全管面具有重要作用.

圖5 管內(nèi)流場(chǎng)溫度等勢(shì)分布圖（單位：K）

圖6為流場(chǎng)壓力等勢(shì)分布圖，超音速?gòu)椡栾w行過程中，首先在前楔表面形成激波壓縮區(qū)，導(dǎo)致前楔表面的壓力較高（p＝1 MPa），對(duì)彈丸的飛行產(chǎn)生極大阻力.前楔斜激波在管壁發(fā)生反射，氣體被再次壓縮，壓力急劇升高，在管道喉部附近的管壁處產(chǎn)生高壓區(qū)（壓力峰值為2.07 MPa）.但是，氣流通過管道喉部時(shí)，截面迅速增大，形成典型普朗特－邁耶膨脹，結(jié)束前楔反射激波波后高壓區(qū)，并在彈丸肩部后方形成低壓區(qū)（p≤0.15 MPa），對(duì)彈丸的飛行產(chǎn)生極大阻力.而在后楔表面的反射作用下，后楔表面壓力恢復(fù)至前楔表面的壓力水平，所以，前楔反射激波在后楔表面的反射位置，影響彈丸表面的壓力分布，并進(jìn)一步影響彈丸的推力（或阻力），應(yīng)盡量控制反射點(diǎn)靠近彈丸肩部，以增加對(duì)彈丸的推力（或減小阻力）.在彈丸底部，氣流通過的截面增大，發(fā)生普朗特－邁耶膨脹，使得底部壓力相對(duì)較低，同時(shí)，由于復(fù)雜波系與旋渦的作用，彈丸底部上下兩側(cè)的壓力并不均衡，這可能影響彈丸飛行的穩(wěn)定性及系統(tǒng)的整體性能.

圖6 管內(nèi)流場(chǎng)壓力等勢(shì)分布圖（單位：MPa）

3 結(jié)論

基于二維非定?？蓧篍uler方程，通過高精度WENO格式與自適應(yīng)加密笛卡爾網(wǎng)格，數(shù)值模擬了沖壓加速器冷態(tài)實(shí)驗(yàn).計(jì)算結(jié)果與實(shí)驗(yàn)[11]一致性很好，證明了數(shù)值方法的有效性.計(jì)算結(jié)果清晰地展示了流場(chǎng)中波系的結(jié)構(gòu)和渦的生成與發(fā)展、激波與渦街的相互作用，最后形成波系結(jié)構(gòu)復(fù)雜的管內(nèi)流場(chǎng)的整個(gè)過程，同時(shí)，給出無粘條件下，彈丸周圍及下游一段距離范圍內(nèi)的重要物理參數(shù)——溫度與壓力，而激波結(jié)構(gòu)、溫度與壓力的分析對(duì)優(yōu)化彈丸結(jié)構(gòu)、速度、材料分布和加速管內(nèi)裝填壓力有重要意義，可以給沖壓加速器的冷態(tài)實(shí)驗(yàn)提供有益的指導(dǎo).

[1]HERTZBERG A，BRUCKNER A P，BOGDANOFF D W.Ram accelerator：a new chemical method for accelerating projectiles to ultrahigh velocities[J].AIAA，1988，26（2）：195－203.

[2]BRUCKNER A P，BOGDANOFF D W，KNOWLEN C，et al.Investigation of gasdynamic phenomena associated with the ram accelerator concept，AIAA－87－1327[R].1987.

[3]HERTZBERG A，BRUCKNER A P，KNOWLEN C.Experimental investigation of ram accelerator propulsion modes[J].Shock Waves，1991，1：17－25.

[4]BENGHERBIA T，YAO Y，BAUER P，et al.Numerical investigation of thermally choked ram accelerator in sub－detonative Regime，AIAA 2009－635[R].2009.

[5]BENGHERBIA T，YAO Y，BAUER P，et al.Thrust prediction in thermally choked ram accelerator，AIAA2010－1129[R].2010.

[6]柳森，白智勇，簡(jiǎn)和祥.沖壓加速器的冷發(fā)射試驗(yàn)[J].流體力學(xué)實(shí)驗(yàn)與測(cè)量，1997，11（4）：8－12.LIU Sen，BAI Zhi－yong，JIAN He－xiang.Ram accelerator cold shot of RAMAC37[J].Experiments and Measurements in Fluid Mechanics，1997，11（4）：8－12.（in Chinese）

[7]柳森，簡(jiǎn)和祥，白智勇，等.37mm沖壓加速器試驗(yàn)和計(jì)算[J].力學(xué)學(xué)報(bào)，1999，31（4）：450－455.LIU Sen，JIAN He－xiang，BAI Zhi－yong，et al.Experimental tests and numerical calculations for the 37mm ram accelerator[J].Acta Mechanica Sinica，1999，31（4）：450－455.（in Chinese）

[8]陳堅(jiān)強(qiáng)，張涵信，高樹椿.沖壓加速器燃燒流場(chǎng)的數(shù)值模擬[J].空氣動(dòng)力學(xué)學(xué)報(bào)，1998，16（3）：297－303.CHEN Jian－qiang，ZHANG Han－xin，GAO Shu－chu.Numerical simulation of the supersonic combustion flowfield in a ram accelerator[J].Acta Aerodynamica Sinica，1998，16（3）：297－303.（in Chinese）

[9]張國(guó)強(qiáng)，金志明，翁春生.沖壓加速機(jī)理數(shù)值研究[J].彈道學(xué)報(bào)，2000，12（3）：27－31.ZHANG Guo－qiang，JIN Zhi－ming，WENG Chun－sheng.Numerical study of principles of ram accelerator[J].Journal of Ballistics，2000，12（3）：27－31.（in Chinese）

[10]張國(guó)強(qiáng)，金志明，翁春生.沖壓加速器燃燒流場(chǎng)數(shù)值模擬[J].彈道學(xué)報(bào)，2001，13（3）：38－41.ZHANG Guo－qiang，JIN Zhi－ming，WENG Chun－sheng.Numerical simulation of combustion flow－field in ram accelerator[J].Journal of Ballistics，2001，13（3）：38－41.（in Chinese）

[11]YATSUFUSA T，CHANG X，TAKI S.Experiments on flame holding position of the finless projectile in ram accelerator，AIAA 2001－1765[R].2001.

[12]LIU Xu－dong，OSHER S，CHAN T.Weighted essentially nonoscillatory schemes[J].J Comput Phys，1994，115：200－212.

[13]BERGER M，COLELLA P.Local adaptive mesh refinement for shock hydrodynamics[J].J Comput Phys，1988，82：64－84.

[14]MITTAL R.Immersed boundary method[M].Annu Rev Fluid Mech，2005，37：239－261.

[15]SCHARDIN H.High frequency cinematography in the shock tube[J].J Photo Sci，1957，5：19－26.