基于Matlab/Simulink的超聲速自由射流試驗系統建模分析

孫順利 李綱 蘆海洋

摘 要： 為實現對關鍵指標艙壓的動態變化過程仿真， 基于Matlab/Simulink對某超聲速自由射流試驗系統的關鍵物理過程建模， 并將仿真結果與試驗結果比較分析。該模型的仿真艙壓與試驗艙壓的曲線變化規律一致， 試驗狀態穩定段的仿真艙壓與試驗艙壓的誤差在1 kPa以內， 噴管出口仿真靜壓值與試驗靜壓值也小于1 kPa， 同時還具備模擬啟動遲滯現象和啟動后的艙壓與噴管進口總壓的正線性關系現象的能力。對比分析結果表明， 該Simulink建模方法能有效模擬超聲速自由射流試驗系統艙壓的動態變化。

關鍵詞：自由射流試驗;? 艙壓;? 超聲速噴管;? 引射;? 擴壓

中圖分類號：??? TJ763; V231.3 ?文獻標識碼：??? A? 文章編號：1673-5048（2021）04-0076-06

0 引? 言

超聲速自由射流試驗系統可在地面模擬高空飛行條件， 常用于沖壓發動機的考核試驗， 是研制沖壓發動機的基礎設施之一， 吸引了大量科研機構進行研究[1-4]。其工作原理是加熱后的高壓氣體通過自由射流噴管達到超聲速狀態， 對高空艙內的氣體產生引射抽吸作用， 使高空艙內的壓力下降， 同時自由射流噴管出口氣流充分膨脹， 獲得滿足沖壓發動機進氣道進口前的均勻馬赫數和壓力分布， 發動機試驗件放置在該均勻流場內。超聲速氣流與沖壓發動機的排氣混合后， 在下游擴壓器減速增壓作用和主動引射器氣流抽吸作用下， 經亞聲速擴壓器排出到大氣環境。 關鍵部件和結構關系如圖1所示。

對于超聲速自由射流試驗系統， 高空艙模擬壓力（簡稱艙壓）指標直接反映自由射流噴管的啟動狀態， 是試驗人員判斷試驗系統工作狀態的重要指標。文獻[5]采用一維理論方法計算了某型超聲速試驗系統的性能， 指出艙壓最小值由擴壓器壅塞和射流膨脹現象決定， 且試驗件堵塞會導致艙壓升高;? 文獻[6]采用試驗方法研究了真空引射的啟動與不啟動現象的運行機理和影響因素;? 文獻[7-8]綜合運用理論分析、 試驗和數值仿真方法深入研究了超聲速試驗系統的啟動特性、 負載匹配和壓力恢復性能等， 獲得設備結構和氣動參數對艙壓的影響;? 文獻[9-10]研究了超聲速引射器啟動遲滯現象的影響因素和預測方法。上述文獻大都側重于單個穩定狀態點的性能指標分析， 對艙壓指標的連續動態變化也多采用試驗方法， 而采用理論方法進行研究的較少。本文借鑒文獻[11]提出的作用力平衡的研究方法和文獻[12]中的建模方法， 基于Matlab/Simulink對某超聲速自由射流試驗系統的關鍵物理過程建模， 實現對試驗系統艙壓的連續動態變化過程仿真， 并與試驗結果比較分析， 驗證建模方法的可行性。

1 模型建立

1.1 高空艙數學模型

根據質量守恒定律， 高空艙控制體流入流出的質量流量總變化等于高空艙內的質量變化， 進而可計算出高空艙內的質量， 再通過理想氣體方程計算出艙壓：

mc=mcini+∫（m·leak- m·s）dt=PaRTVc+∫（m·leak- m·s）dt（1）

Pc=mcVcRT （2）

高空艙流入流量主要是高空艙補氣造成， 高空艙補氣過程簡化為限流孔節流絕熱過程， 當艙內和艙外的壓差越大， 則補氣量越大， 可得補氣流量[12]： 航空兵器 2021年第28卷第4期

孫順利， 等： 基于Matlab/Simulink的超聲速自由射流試驗系統建模分析

m·leak=

CfAPa2γRT（γ-1）PcPa2γ-PcPa）γ+1γ?? PcPa≥2γ+1γ+1γ-1

CfAPaγRT2γ+1γ+1γ-1

PcPa<2γ+1γ+1γ-1（3）

式中： mc和mcini為高空艙內氣體的實時和初始質量;? m·leak和m·s為高空艙補氣和流出的質量流量;? Pa和Pc為大氣壓和高空艙實時靜壓;? R， T和γ為氣體常數、 靜溫和比熱比;? Vc為高空艙容積;? Cf和A為補氣孔流量系數和面積。? 其中補氣孔面積和流量系數可根據試驗和調試數據進行修正。

1.2 自由射流引射數學模型

高空艙流出流量由自由射流引射數學模型決定，射流噴管出口的靜壓P1n大于噴管出口艙壓P1s（艙壓等于二次流靜壓）， 即當P1n>P1s， 噴管出口氣流等熵膨脹：

P2n=P1nν1nν2nγ=P1nAfnAshd-As-Ablkγ（4）

二次流氣流等熵壓縮：

P2s=P1s1+γ-12M22s1+γ-12M21sγγ-1（5）

兩股氣流在下游某點達到壓力平衡（P2n=P2s）， 形成氣動喉道As：

As=Ashd-Ablk-AfnP1sP1n1+γ-12Ma22s1+γ-12Ma21sγγ-11γ （6）

式中： Ablk為試驗件等堵塞面積;? γ為比容;? 當P1n≤P1s時， 噴管出口氣流不膨脹， 氣動喉道為擴壓器管道面積減去自由射流噴管出口面積， 即As=Ashd-Afn， 此處不包含發動機試驗件和臺架， 因而不存在堵塞面積。

（1）當As處達到壅塞時， Ma2s=1， 可得二次流的最大質量流量m·smax：

m·smax=γR2γ+1γ+1γ-1PcTcAs （7）

通過As與A1s的等熵膨脹關系， 可得二次流進口的馬赫數Ma1s和靜壓P1s， 最后得到作用力Fsmax：

Fsmax=P1sA1s（1+γMa21s）（8）

需要注意的是， 當As≥0時， m·smax≥0， 此時二次流流出高空艙控制體外;? 而當As<0時， m·smax<0， 射流噴管氣流在擴壓器內欠膨脹， 此時二次流流入高空艙控制體內， 艙壓升高， 這是模擬當自由射流噴管啟動后， 艙壓隨自由射流噴管進口總壓升高而升高的線性關系現象[13-14]的關鍵方法。

（2）當As處未壅塞時， 假設二次流流量在正負最大二次流流量之間， 即-m·smax

Ma1s-m·sA1sPtγRTt1+γ-12Ma21sγ+12（γ-1）=0（9）

采用二分法迭代求解， 得亞聲速解Ma1s（0～1）。 由等熵膨脹關系得靜壓P1s：

P1s=Pc1+γ-12Ma21s-γγ-1? （10）

最后得二次流作用力Fs：

Fs=P1s（Ashd-Afjn）（1+γMa21s）（11）

式中： 1， 2， mix為擴壓器不同截面;? n， s為自由射流和二次流;? Ma， Pt， Tt為馬赫數、 總壓和總溫。

1.3 自由射流噴管氣動數學模型

拉瓦爾噴管出口氣流靜壓和馬赫數受進口總壓和出口背壓的共同影響， 存在9種情況[15]。 根據仿真建模需要可簡化為以下3種情況：

（1）當π1

P1n=Pc（12）

Ma1n=2γ-1（1-（P1n/Pt）γ-1-γ）（13）

（2）當π2≤Pc/Pt≤π1， 喉道壅塞， 噴管擴張段內存在正激波， 噴管出口靜壓為艙壓， 由喉道壅塞條件和質量守恒定律可得

P1n=Pc（14）

Ma1n=1γ-1+

1γ-12+2γ-12γ+1γ+1γ-1PtP1n2AtA1n2（15）

式中： At和A1n為噴管喉道和出口面積。? 出口背壓與進口總壓的臨界壓比π1與π2的計算方法見文獻[16-17]。

（3）當Pc/Pt<π2， 喉道壅塞， 擴張段內全部為超聲速， 出口馬赫數為噴管名義馬赫數Mafn， 由噴管喉道和出口面積比決定。 由等熵膨脹關系式可得出口靜壓：

Ma1n=Mafn（16）

P1n=Pt（1+γ-12Ma2fn）-γγ-1（17）

最后將噴管出口馬赫數Ma1n和靜壓P1n帶入下式， 得到質量流量和作用力：

m·n=Ma1nA1nPtγRTt/（1+γ-12Ma21n）γ+12（γ-1）（18）

Fn=P1nA1n（1+γMa21n）（19）

1.4 超聲速擴壓器擴壓過程數學模型

擴壓器出口壓力主要基于質量、 能量和動量守恒對自由射流和二次流的氣流混合求解出口參數：

m·n+m·s=m·mix（20）

m·nCpnTtn+m·sCpnTts=m·mixCpmixTtmix（21）

P1nA1n+P1sA1s-PmixAmix-Ff=m·mixvmix-m·nvn-m·svs（22）

式中： Cp， v為定壓比熱容和速度。

假設管道摩擦力Ff=0， 作用力的計算式為

F=PA+m·v=PA+PRTAv2=PA（1+γγRTv2）=

PA（1+γMa2）（23）

替換式（22）中相關項后可得

F1n+F1s=Fmix（24）

最后可得到擴壓器出口的質量流量、 總溫和作用力。

質量流量計算式為

m·=PAMaγRTt（1+γ-12M2）（25）

式（25）與式（23）中靜壓P相等， 聯立后得

P=m·AMγRTt（1+γ-12M2）=FA（1+γM2） （26）

經整理得關于M2二次方程式：

AM4+BM2+C=0（27）

式中： A，B，C分別表示代替二次方程系數。

由于混合過程具備真實物理意義， 因此B2-4AC≥0， 式（27）在數學意義上必然存在實數解：

M2=-B±B2-4AC2A（28）

式中： 取+號時為超聲速解， 當該值為負數， 則不存在超聲速解;? 取-號時為亞聲速解， 對應于存在超聲速解時的正激波波后亞聲速解或者不存在超聲速解時的亞聲速解。 當特殊情況B2-4AC=0， 此時恰好M=1。

最后可得混合后的擴壓器出口靜壓Pout， 考慮到由于在超擴段內通常不是理想正激波， 而是以多道斜激波的激波串使壓力恢復， 因此存在壓力損失， 考慮超聲速擴壓器內非理想正激波壓力損失修正的的波后靜壓：

Pmix-crt=Pmix（1-σγmix-12Ma2mix）γmixγmix-1 （29）

考慮亞聲速擴壓器內的膨脹作用的出口壓力：

Pout=Pmix-crt（1+ηdγmix-12Ma2mix）γmixγmix-1（30）

式中： σ為超聲速段內的壓力恢復系數;? ηd為擴壓效率。

1.5 動態仿真的原理和基本流程

高空艙流出流量由自由射流引射數學模型決定， 對應于二次流m·s， 其數學模型建立過程如圖2所示。 動態仿真模型流程如圖3所示。其關鍵算法是計算出正確的二次流流量， 當混合后的擴壓器出口壓力大于或等于大氣壓時， 最大的二次流流量即為二次流流量， 當混合后的出口壓力小于大氣壓時， 則二次流量介于正最大二次流流量和負最大二次流流量之間， 通過二分法求得實際二次流流量， 使出口壓力等于大氣壓。

圖4是H=6 km， Ma=2.5工況下， 二次流最大值m·smax和實際值m·s的對比。當m·s

2 仿真與試驗結果對比分析

圖5～6是自由射流試驗系統的兩次試驗及仿真結果。需要說明的是， 此高度范圍不需要打開下游的引射器， 實際模擬高度以試驗結果的靜壓為準。噴管出口靜壓的測量點位于噴管出口內壁面。

對比圖5～6， 可見仿真艙壓和試驗艙壓的變化規律基本一致， 對噴管進口總壓參數的變化具有較好的響應跟隨性， 當自由射流噴管總壓和艙壓穩定后， 仿真艙壓與試驗艙壓的誤差在1? kPa以內，? 該Simulink模型可以準確模擬試驗系統穩定后艙壓值。

另外， 兩次試驗仿真中， 當自由射流噴管總壓到目標值時（圖5（b）中17～21 s， 圖6（b）中13～50 s）， 補氣孔面積差異導致艙壓高于或低于相應模擬高度， 當補氣孔面積足夠大時， 艙壓無法抽到模擬高度附近;? 補氣孔面積足夠小時， 艙壓遠小于模擬高度。因此， 補氣孔面積是艙壓調節的關鍵參數。

當艙壓受補氣孔面積影響而沒有與實際模擬高度一致時， 自由射流噴管出口靜壓的試驗值和仿真值均達到模擬高度， 二者誤差在1 kPa以內， 噴管出口達到名義馬赫數為2.5。 結果表明在艙壓高于和低于噴管出口靜壓的兩種條件下， 該Simulink模型也可以準確模擬試驗系統穩定后的噴管出口靜壓條件。因此， 仿真模型通過準確計算的噴管出口艙壓和靜壓， 根據斜激波理論和普朗特邁耶理論可計算噴管出口的均勻流場的菱形區面積[18]， 獲得發動機試驗件的安裝位置和有效試驗時間。

艙壓與自由射流噴管進口總壓變化關系對比如圖7所示。在啟動過程中， 艙壓先隨自由射流噴管進口總壓降低， 當艙壓達到最低時， 隨自由射流噴管總壓呈線性關系增加， 與文獻[14]中描述一致。在圖7（a）中， 仿真結果也顯示艙壓會隨自由射流噴管總壓變化時出現啟動遲滯現象。 與試驗結果相比， 圖7（a）中的仿真最小啟動壓力（A點）與試驗最小啟動壓力（B點）相差在70 kPa左右， 仿真和試驗的最小保持啟動壓力（C點）基本吻合， 可見仿真結果較好地捕捉到啟動遲滯現象的兩個關鍵試驗點。對于自由射流噴管進口閥門開啟和關閉的動態過程中， 仿真艙壓與試驗艙壓的啟動遲滯幅度存在較大差異， 原因是： （1）自由射流噴管進口閥門的關閉時間比開啟時間要短， 造成噴管進口總壓在下降時的速度比上升時快， 進而在計算二次流流量時， 噴管進口閥門關閉時的實際二次流量較大， 并且容易受到負最大二次流量的限流（m·s=-m·smax）， 而噴管進口閥門開啟時的實際二次流量小， 并且不易造成正最大二次流量的限流（m·s<+m·smax）， 結果可以在圖4中14～16 s和23～24 s時間段觀察到。（2）查閱文獻[19]， 啟動遲滯現象普遍存在亞聲速與超音速轉變過程， 如進氣道和引射器的啟動與不啟動狀態轉變過程。 啟動時由亞聲速轉變為超聲速需要突破正激波的壓力損失， 因而需要更高的進口總壓， 而從超聲速轉變為亞聲速時不啟動狀態則壓力損失較小， 特別是對于存在二次喉道的情況， 這種啟動遲滯現象會更明顯。當前仿真模型中超聲速擴壓段的壓力恢復系數σ是固定不變的， 圖7（b）是通過在啟動和關閉階段設置不同的壓力恢復系數的仿真結果， 可見在噴管進口閥門開啟和關閉階段， 曲線的啟動遲滯差異縮小。結果表明超聲速壓力恢復系數會對仿真艙壓結果產生較大影響， 因此在實際應用中需要根據試驗結果對不同試驗時間段內的壓力恢復系數σ修正， 合理取值。需要注意的是， 在圖5～6中， 當仿真艙壓與自由射流噴管進口總壓超過臨界壓比π2=0.417時， 自由射流噴管出口從亞聲速（Ma=0.51）突變為超聲速（Ma=2.5）， 即將噴管擴張段內的正激波從出口推了出去， 自由射流噴管達到滿流狀態， 同時噴管出口的仿真靜壓也突然下降至最低， 然后噴管出口的仿真靜壓隨著自由射流噴管總壓的升高而升高， 最終穩定在可供試驗吹風的狀態， 表明噴管數學模型較好地模擬了噴管理想狀態的啟動過程。另外， 注意到在噴管啟動過程中， 噴管出口仿真靜壓與試驗靜壓出現差異， 主要由于試驗靜壓測量位置位于靠近噴管出口內壁面處， 該位置在艙壓遠高于噴管出口靜壓時， 極易出現邊界層分離， 而當前理想噴管數學模型假設噴管出口為均勻氣流， 未能考慮噴管出口的邊界層分離的不均勻出口參數， 后續工作需增加噴管模型邊界層分離模擬能力。

3 結? 論

（1） 仿真艙壓與試驗艙壓的曲線變化規律一致， Simulink模型準確地模擬了整個試驗過程中艙壓隨自由射流噴管進口總壓的變化過程， 還準確模擬了自由射流噴管的啟動過程、 艙壓的啟動遲滯現象和艙壓與自由射流噴管進口總壓的正線性關系等現象， 并且試驗狀態穩定段的噴管出口仿真靜壓與試驗靜壓高度吻合， 說明應用Simulink建模方法模擬超聲速自由射流試驗系統性能具備可行性和工程應用性。

（2） 補氣孔面積是艙壓調節的關鍵參數， 試驗前需要準確預估補氣閥參數設置。壓力恢復系數對啟動遲滯現象影響較大， 需要根據試驗結果對不同試驗時間段的壓力恢復系數合理修正， 提高仿真準確度。

（3） 目前射流噴管數學模型未包含邊界層分離現象， 在噴管啟動過程中的噴管出口的仿真值與試驗值存在差異， 也沒有將引射器數學模型和發動機排氣數學模型考慮在內， 在后續工作中， 如高空高馬赫數建模分析時， 將會增加這部分內容， 提高仿真準確度。

參考文獻：

[1] Dunsworth L C， Reed G J. Ramjet Engine Testing and Simulation Techniques[J]. Journal of Spacecraft and Rockets， 1979， 16（6）：? 382-388.

[2] Marren D， Lu F. Advanced Hypersonic Test Facilities[M]. Reston， VA：? AIAA， 2002.

[3] 樂嘉陵. 吸氣式高超聲速技術研究進展[J]. 推進技術， 2010， 31（6）：? 641-649.

Le Jialing. Progress in Air-Breathing Hypersonic Technology[J]. Journal of Propulsion Technology， 2010， 31（6）：? 641-649.（in Chinese）

[4] 韓建濤， 孫順利， 李綱， 等. 固沖發動機進氣道半自由射流試驗直管擴壓器研究[J]. 航空兵器， 2019， 26（4）：? 95-98.

Han Jiantao， Sun Shunli， Li Gang， et al. Straight Divergent Tube Research of Semi-Freejet Experiment for Ramjet Inlet[J]. Aero Weaponry， 2019， 26（4）：? 95-98.（in Chinese）

[5] Nagaraja K S， Hammond D L， Graetch J E. One-Dimensional Analysis of Compressible Ejector Flows Applicable to V/STOL Aircraft Design[EB/OL].（2012-08-16）[2020-09-02]. AIAA Journal.https：∥doi.org/10.2514/6.1973-1184.

[6] Arun K R， Rajesh G. Flow Transients in Un-Started and Started Modes of Vacuum Ejector Operation[J]. Physics of Fluids， 2016， 28（5）：? 056105.

[7] 吳繼平. 高增壓比多噴管超聲速引射器設計理論、 方法與實驗研究[D]. 長沙：? 國防科學技術大學， 2007.

Wu Jiping. Design Theory， Method and Experimental Investigation of High Compression Ratio Multi-Nozzle Supersonic Ejector[D]. Changsha：? National University of Defense Technology， 2007. （in Chinese）

[8] 陳健. 超-超引射器內部流動過程研究[D]. 長沙：? 國防科學技術大學， 2012.

Chen Jian. Researches on the Flow Process of the Supersonic-Supersonic Ejector[D]. Changsha：? National University of Defense Technology， 2012. （in Chinese）

[9] Park G， Kim S， Kwon S. A Starting Procedure of Supersonic Ejector to Minimize Primary Pressure Load[J]. Journal of Propulsion and Power， 2008， 24（3）：? 631-635.

[10] Kim S， Kwon S. Starting Pressure and Hysteresis Behavior of an Annular Injection Supersonic Ejector[J]. AIAA Journal， 2008， 46（5）：? 1039-1044.

[11] Daniel D. A General Simulation of an Air Ejector Diffuser System[C]∥28th Aerodynamic Measurement Technology， Ground Testing， and Flight Testing Conference， 2012.

[12] 馬飛， 王海洲. 基于Matlab/Simulink的氣體增壓系統建模分析[J]. 導彈與航天運載技術， 2006（3）：? 41-47.

Ma Fei， Wang Haizhou. Modeling Analysis of a Gas Pressurized System Based on Matlab/Simulink[J]. Missiles and Space Vehicles， 2006（3）：? 41-47.（in Chinese）

[13] Kim S， Kwon S. Experimental Investigation of an Annular Injection Supersonic Ejector[J]. AIAA Journal， 2006， 44（8）：? 1905-1908.

[14] 孫順利， 李綱. 高空模擬引射器啟動特性的實驗和數值計算研究[J]. 彈箭與制導學報， 2017， 37（4）：? 63-67.

Sun Shunli， Li Gang. Experimental and Numerical Research on? the Starting Characteristics of the High Altitude Simulation Ejector[J]. Journal of Projectiles，Rockets，Missiles and Guidance， 2017， 37（4）： ?63-67.（in Chinese）

[15] White F M. Fluid Mechanics[M]. New York：? McGraw-Hill， 2003： 599-632.

[16] Anderson J D. Modern Compressible Flow：? With Historical Perspective[M]. New York：? McGraw-Hill， 1990.

[17] 周文祥， 黃金泉， 周人治. 拉瓦爾噴管計算模型的改進及其整機仿真驗證[J]. 航空動力學報， 2009， 24（11）：? 2601-2606.

Zhou Wenxiang， Huang Jinquan， Zhou Renzhi. Improvement of Laval Nozzle Calculation Model and Simulative Verification in Aero-Engine Performance Calculation[J]. Journal of Aerospace Power， 2009， 24（11）：? 2601-2606.（in Chinese）

[18] Pruitt D， Bates L. Starting and Test Rhombus Characteristics of Two-Dimensional Supersonic Free-Jet Nozzle/Generic Supersonic Aircraft Inlet Configurations[C]∥AlAA 4th International Aerospace Planes Conference， 1992.

[19] Van Wie D M， Kwok F T， Walsh R F. Starting Characteristics of Supersonic Inlet[R]. AIAA 1996-2914.

Modeling Analysis of Supersonic Free Jet Test System

Based on Matlab/Simulink

Sun Shunli ， Li Gang， Lu Haiyang

（China Airborne Missile Academy， Luoyang 471009， China）

Abstract：?? In order to simulate the dynamic change process of the key index cabin pressure， the key physical process of supersonic free jet test system is modeled based on Matlab/Simulink， and the simulation results are compared with the test results. The curve variation law of the simulated cabin pressure of the model is consistent with that of the test cabin pressure. The error between the simulated cabin pressure and the test cabin pressure in the stable section of the test state is within 1 kPa， and the simulated static pressure at the nozzle outlet and the test static pressure are also less than 1 kPa. At the same time， it also has the ability to simulate the start-up hysteresis phenomenon and the main line relationship between the cabin pressure after startup and the total pressure at the nozzle inlet. The comparative analysis results show that the Simulink modeling method can effectively simulate the dynamic change of cabin pressure of supersonic free jet test system.

Key words：?? free jet test;? cabin pressure;? supersonic nozzle;? ejector;? diffuser