基于立方型氣體狀態方程的激波風洞準一維流動數值研究

摘要： 采用基于立方型氣體狀態方程的準一維流動數值模擬方法研究了反射式高焓激波風洞的真實氣體流動，重點關注了高壓真實氣體效應對風洞全場流動時空結構和駐室區氣流參數的影響，并以理論分析揭示了高壓真實氣體效應對激波管內流動的作用機理。研究表明：對于以冷高壓氣體驅動的激波風洞，使用考慮分子體積和分子間作用力的真實氣體狀態方程能夠更準確地描述氣體的狀態和風洞內的流動狀況。高壓真實氣體效應主要在冷驅動氣體中發生作用，其作用效果主要是使當地聲速增大，從而使得入射稀疏波和反射稀疏波的傳播速度加快；另一方面，高壓氣體效應在高溫氣體效應較顯著的被驅動氣體中作用微弱，且對激波管產生激波的強度和激波后的流動狀態影響甚微。稀疏波的加快傳播改變了激波管波系的相干時空關系。提前抵達的稀疏波可在一定情況下侵蝕激波風洞的有效試驗時間。對于所測試的激波風洞構型，在150 MPa 氫氣驅動110 kPa 氮氣的工況下，高壓效應導致的有效試驗時間縮短約38%。適當加長驅動段長度和采用高溫氣體驅動均可有效減弱高壓真實氣體效應的影響。

關鍵詞： 激波風洞；高壓真實氣體效應；立方型氣體狀態方程；準一維數值模擬；有效試驗時間

中圖分類號： O354.7 國標學科代碼： 13025 文獻標志碼： A

風洞試驗是認識高超聲速流動機理和預測飛行器性能不可或缺的手段。激波風洞是支撐高馬赫數飛行地面氣動試驗研究的主力設備之一。它主要利用運動激波對流體的氣動加速和氣動壓縮作用產生試驗所需的高焓氣流，因此，其工作效能極大地依賴于噴管上游激波管內以運動激波、稀疏波等為特征的強非定常流動。對激波風洞整體流動的充分理解和認知是此類風洞有效服務于氣動試驗的重要前提。

為了達到模擬高馬赫數飛行環境所需的總溫、總壓條件，激波風洞運行時，驅動段以及激波反射后的駐室區域的氣體往往處于高溫、高壓的狀態。在極端的溫度和壓強條件下，氣體的物理化學屬性與狀態顯著偏離理想氣體的描述，從而在流動中產生所謂的真實氣體效應。這里的真實氣體效應實際包括2 類。一是高溫真實氣體效應。由于氣體溫度的提升，氣體分子原本大量處于基態的振動能和電子能被顯著激發，由此產生分子多種內能模態間（轉動能、平動能、振動能和電子能）有限速率的相互轉化；同時，熱運動加劇使得分子內部化學鍵斷裂乃至電子脫離原子核束縛，氣體發生離解和電離。這種高溫氣體中發生的熱弛豫和化學反應，導致氣體宏觀可感內能的改變以及氣體物理化學屬性的變化[1]。二是高壓真實氣體效應。在氣體溫度維持低值時，高壓意味著分子數密度增加；當氣體狀態逼近液化臨界或超臨界態，氣體分子體積及分子間作用力逐漸變得不可忽略，經典氣動理論使用的理想氣體狀態假設將不再適用，氣體的熱容、內能、熵等熱力學參數同時發生一定程度的偏移[2]。高狀態運行的激波風洞流動中常常兼具上述2 類真實氣體效應，基于經典氣體動力學理論的計算已不足以準確評估風洞的運行狀態。

對激波風洞流動的測試和研究是高超聲速試驗技術的重要組成部分，貫穿于風洞的設計、運行和對試驗結果的分析環節，如：Jacobs[3] 發展了準一維數值模擬程序L1D，用于脈沖風洞內流場模擬，并成功服務于昆士蘭大學T4 風洞運行狀態的預測評估；Boyce 等[4] 對T4 激波風洞不同工況下的試驗氣流參數以及有效試驗時間進行了分析；Hannemann[5] 和Johnston 等[6] 以試驗結合數值模擬考察了HEG 高焓激波風洞復現X-38 等飛行器再入大氣條件的流動狀態；Holden 等[7] 在LENS-I 和LENS-II 激波風洞中獲得了高總溫和高總壓的試驗氣流，并評估了真實氣體對超高速再入飛行器性能的影響；盧洪波等[8] 在FD-21 高焓激波風洞中開展高馬赫數吸氣推進試驗技術探索。近20 年來，人們對于高焓脈沖風洞的流動模擬和計算普遍注意將高溫真實氣體效應納入考慮。Llobet 等[9] 考慮振動非平衡和化學平衡對激波風洞的運行過程進行計算，獲得試驗氣流參數和有效試驗時間隨初始壓比的變化規律。2022 年，Lynch 等[10]使用L1D4 程序（耦合NASA-CEA[11] 計算熱化學平衡）為自由活塞驅動反射型激波風洞設計了總壓20 MPa、總溫3 700 K 的工況，復現了來流馬赫數為8～9 的飛行環境；Luo 等[12] 采用準一維數值模擬方法，考慮化學平衡模型，對JFX 爆轟驅動激波風洞內的流動進行精確快速的模擬和研究，掌握了激波風洞運行過程中的波系運動情況。但學者們對高焓脈沖風洞中高壓真實氣體效應的研究則較少。MacLean 等[13-14] 在對LENS-I 激波風洞和LENS-XX 膨脹管風洞的分析過程中使用了考慮體積修正的狀態方程和熱平衡氣體，發現考慮體積修正后的自由流單位雷諾數明顯增大；Candler[15] 在使用克勞修斯狀態方程和熱化學非平衡的基礎上，對AEDC 激波風洞的噴管流動進行計算分析，并指出克勞修斯狀態方程在低溫高密度狀態下無法完全體現高壓氣體效應。由于考慮完全的立方型狀態方程后，流場溫度、壓強和密度的非線性關系顯著增加了數理建模與流動解算的復雜性，Sirignano[16] 提出對成熟的真實氣體狀態方程進行線性化處理，并用線性化后的狀態方程分析了等熵流等基本流動過程，但是這些簡化處理僅在有限的條件范圍內適用，當流場狀態跨度較大或者測試條件范圍較寬時，簡化方程偏差顯著。因此，仍有必要采用完備的狀態方程，系統地分析激波風洞流動中高壓真實氣體的影響效應，綜合評估考慮高壓氣體效應前后激波風洞波系結構和試驗氣流參數的差異，揭示高壓真實氣體效應的作用機制。

本文中，以FD-14A 激波風洞[17] 為基本構型，基于自主發展的準一維熱化學非平衡數值模擬平臺，引入完全的立方型氣體狀態方程描述高壓真實氣體狀態，對該型風洞在試驗氣流總溫2 000～6 000 K、總壓20～200 MPa 典型工況下的流動過程進行計算和分析；考慮到激波風洞運行過程中高溫真實氣體效應與高壓真實氣體效應是同時存在的，在考慮熱化學非平衡基礎上對比采用真實氣體狀態方程和理想氣體狀態方程的計算結果，著重關注高壓真實氣體效應導致的激波風洞運行全場流動時空特性和試驗氣流參數的差異以及產生差異的物理機制。

1 風洞構型與數值模擬方法

圖1 為FD-14A 激波風洞[17] 內流道的基本構型。該風洞激波管部分管道的內徑均為0.15 m；當前研究設置驅動段長9 m，被驅動段長18 m。為考慮反射端的質量通量，設置噴管喉道的直徑為23.9 mm。

確定喉道和收縮段后，噴管其余部分不會對實際激波管流動產生顯著影響，為了保留激波風洞的全部結構特征，計算過程中保留了整個噴管。

該風洞驅動段氣體為高壓氫氣（H2），被驅動段氣體為空氣或氮氣（N2）。本研究中，為重點考察高壓真實氣體效應，被驅動段統一采用氮氣。

1.1 準一維流動控制方程

對于以大長徑比管道流動為特征的高焓脈沖風洞設備，準一維流動數值模擬一直是一種行之有效且快捷的計算分析手段[3， 12， 14]。通過合理的假設，準一維數值模擬可將管道壁面的摩擦與傳熱、高溫氣體熱化學非平衡等各種物理化學效應同時納入考量。本文中所使用的考慮雙溫度振動非平衡、化學反應和輸運效應的準一維流動控制方程如下：

式中：U 為守恒項；F 為對流項；P 為變截面管道壁面的流向作用力；A 為流道面積函數，描述管道截面空間變化；ρ、u、p、e、ev、Yi 分別為氣體的密度、速度、壓強、質量總能量、振動能和組分 的質量分數。非平衡氣體的質量總能量有：

式中：Ru 為普適氣體常數，wi 為組分i 的分子摩爾質量，a1，i、a2，i、a3，i、a4，i 和a5，i 為該組分的熱力學系數。Δep 為高壓真實氣體效應的能量修正項，與真實氣體狀態方程相關，后續章節將給出推導形式。

方程（1） 中Jin 和Jbd 分別為內流和壁面輸運源項。其中，Jbd 用于描述管道壁面摩擦與傳熱對平均流動的影響。本文中主要采用Groth 等[18] 的壁面效應模型，組合Jacobs 的參考溫度形式，并進行了適配二溫度（平動和振動溫度）非平衡模型的改良。S 為熱化學非平衡源項，包括化學反應組分生成源項和振動能弛豫源項。以上3 個源項的具體形式為：

式中：λi 為由組分擴散導致的第i 組分含量變化，τ 為由黏性引起的流體內部應力，qtr 和qv 分別表示平動能以及振動能的熱傳導，M 表示壁面黏性導致的動量傳遞，Qtr 和Qv 分別表示壁面傳熱導致的平動能和振動能傳遞， 表示由化學反應導致的第i 組分生成，sv 表示由振動弛豫導致的能量轉化。對于計算過程中需使用的氣體輸運系數（主要是黏性系數和熱傳導系數），單一組分的輸運系數以Sutherland 公式計算，混合氣體輸運系數由單一組分系數按Wilke[19] 的方法混合得到。

壁面動量損失源項M 和壁面能量損失源項Q 的表達式分別為：

式中：D 為管道水力學直徑；f 為Darcy-Weisbach 摩擦因數[18]；cp 為流體的比定壓熱容；Pr 為普朗特數；Tw 為當地壁溫；Taw 對平動能為絕熱壁溫，取值為經過可壓縮修正之后的氣流平均溫度，對振動能為當地平均振動溫度。

本文中振動弛豫采用Park[20] 的雙溫度模型，特征弛豫時間由Millikan 等[21] 的半經驗表達式計算得到，并參考Park[20] 的方法對極高溫度下的弛豫時間進行修正。空氣在高溫下的化學反應采用Gupta 等[22]的基元反應動力學機理（本文中僅涉及氫氣和氮氣，因而僅取其中氮氣離解反應）。具體模型和數據可參考文獻[22-23]。

為了在統一模型架構下求解激波風洞的高溫高壓真實氣體流動，在考慮高溫氣體效應的基礎上，采用以下立方型氣體狀態方程：

封閉方程組（1）。式（6） 中：v=1/ρ 為氣體比容，R=Ru/w 為氣體常數，w 為分子摩爾質量；a、b、k1 和k2 為模化分子間相互作用以及分子體積的參數或函數。當a 和b 取值為零時，方程（6） 回歸為理想氣體狀態方程。

1.2 數值計算方法

準一維數值模擬采用基于非結構網格的有限體積法離散并解算方程（1），考慮到源項S 中化學反應的加入會使得方程具有一定的剛性，因此，采用時間分裂步法將流動與熱化學源項解耦求解，首先計算熱化學凍結條件下的流動過程，再在固定網格內求解定容絕熱熱化學弛豫的剛性常微分方程組。純流動過程采用二階時空精度的MUSCL-Hancock 格式求解，數值通量通過HLLC[24] 格式計算得到，時間步長由CFL （Courant-Friedrichs-Lewy） 條件控制；熱化學非平衡源項部分采用成熟的剛性常微分方程求解器DVODE 進行計算。

1.3 數值方法驗證

本研究基于自主開發的準一維熱化學非平衡數值模擬程序，所采用的基本模型和方法已在1.1 節和1.2 節中描述。該程序的基本算法和輸運模型已在多個研究中獲得充分驗證[25-27]，并在實際脈沖風洞結構和工況設計中獲得應用[28]。

首先，采用激波管試驗對上述壁面摩擦和傳熱模型進行校驗。將含壁面輸運模型的準一維數值程序計算得到的激波管中入射激波馬赫數Mas 沿程衰減情況與試驗結果進行對比，如圖2 所示。

圖2 顯示，準一維數值模擬可很好地刻畫廣泛條件下運動激波由于壁面摩擦和傳熱而發生衰減的現象。

接下來，對引入立方型狀態方程后的程序進行補充驗證。選取表1 所示的典型工況，計算激波風洞駐室壓強隨時間演化情況，并與試驗測量壓強數據進行對比，如圖3 所示。

由圖3 可知，數值模擬計算得到激波到達時刻以及駐室區壓強變化情況與試驗數據基本吻合。考慮立方型氣體狀態方程后的數值模擬更準確地捕捉到了稀疏波系的運動情況，計算結果與試驗結果的吻合度顯著提升。以上結果表明，本文中采用的準一維數值模擬方法和程序可以準確預測激波風洞的整體運行狀態。

2 結果與討論

2.1 立方型氣體狀態方程和熱力學參數

狀態方程建立氣體壓強、溫度和密度（或比容、數密度等）之間的關系。理想氣體（ideal gas， IG）狀態方程將氣體分子視為質點，而不考慮氣體分子體積以及分子間作用力，因此在低溫高壓條件下，它無法準確描述氣體狀態。van der Waals[2] 建立了第一個適用于真實氣體的立方型狀態方程，即van derWaals（VDW）方程，該方程將分子體積和分子間作用力的影響分別以斥力項pr 和引力項pa 的形式體現在壓強表達式中[29]：

式中：參數b 表征分子體積，參數a 度量分子間吸引力；考慮到分子體積會對分子間作用力產生影響，因此引力項中還包含關于體積（比容v）的函數g（v），VDW 方程中取g（v）=v2。參數a、b 與氣體臨界參數（臨界溫度Tcr、臨界壓強pcr）相關，其中，a0 和b0 為狀態方程設定的常系數：

后續研究者在VDW 方程基礎上發展出數百種立方型狀態方程[30]，其中，Soave-Redlich-Kwong（SRK） 方程[31] 和Peng-Robinson （PR） 方程[29] 是2 種代表性的立方型氣體狀態方程。他們通過采用形式更復雜的g（v） 并引入溫度修正系數α 來更好地描述氣體在高壓條件下的狀態。

SRK 方程[31] 的具體形式為：

式中：fω 為分子偏心因子ω 的函數。

SRK 方程[31] 和PR 方程[29] 適用于單一組分，對于混合氣體，由于所含各組分的參數a、b 不同，需要在考慮各組分摩爾分數的基礎上進行平均化處理。常見平均化處理方式有以下2 種。一種是基于組分二元交互作用的經典混合規則[32]：

式中：xi、ai 和bi 分別為第i 組分的摩爾分數以及對應的參數，δij 為二元作用參數（這里參考Poling 等[33]提出的最簡單的處理方法，即取δij=0）。Peng 等[29] 在提出PR 方程時對于混合氣體參數采用了上述規則[32]。另一種混合氣體平均化處理方法是Kay[34] 提出的準臨界參數混合規則，該規則首先對混合氣體的臨界溫度和臨界壓強進行平均處理，再由式（8） 計算得到相應的平均參數a 和b。混合氣體的臨界溫度和臨界壓強的計算公式分別為：

式中：Tcr，i 和pcr，i 分別為i 組分的臨界溫度和臨界壓強。

將2 種混合規則計算得到的70% 氫氣和30% 氮氣混合氣體的可壓縮性因子（z=pv/RT）隨氣體壓強p 的變化與NIST 試驗數據[35] 進行對比，如圖4 所示。可以看到，使用經典混合規則計算方法獲得的可壓縮性因子大幅偏離試驗值，而使用準臨界參數混合規則可以較好地描述混合氣體狀態。因此，后續計算將使用式（14） 進行平均處理。

氣體的熱力學參數（內能、焓、熵和比熱等）的具體形式與氣體狀態方程密切相關，氣體狀態方程變化后，需要對上述參數進行重新推導。針對立方型方程一般形式（6），由熱力學基本方程推導獲得真實氣體狀態方程與理想氣體狀態方程內能之差（式（2） 中Δep 項），表達如下：

式中：Δe0，ref=Δe0（v=vref， T=Tref），是與滿足理想氣體假設的參考體積vref 和參考溫度Tref 有關的常數。獲得內能增量后，可進一步計算得到比定容熱容cV、比定壓熱容cp 以及聲速c：

式中：cV，IG 為理想氣體狀態下氣體的比定容熱容，γ 為比熱比。上述方程中保留部分導數以保持方程形式簡潔，這些導數均可由狀態方程本身或 表達式直接求導得出，這里不再展開。

氣體狀態方程的選取，特別是溫度修正系數的具體形式會對氣體的熱力學參數產生不可忽視的影響，選取合適的狀態方程是分析激波風洞時空特性和試驗氣流參數的前提。依據以上各式計算不同狀態方程（VDW 方程、SRK 方程和PR 方程）下氫氣的可壓縮性因子和聲速隨壓強的變化，并與NIST 試驗數據[35] 進行對比，如圖5～6 所示。

3 種立方型方程均能很好地呈現氣體可壓縮性和聲速隨壓強升高逐步升高的特征。其中，VDW 方程在高壓下顯著高估可壓縮性和聲速，而SRK 和PR 方程兩者的結果比較接近，且與NIST 試驗數據較接近。對比后兩者，SRK 方程在可壓縮性預測上更有優勢，而PR 方程對聲速的預測適用范圍更寬。對于壓強150 MPa 的氫氣，PR 方程計算得到的聲速為2 600 m/s，相比試驗值2 480 m/s 高約5%；而此時，理想氣體方程計算得到的聲速為1 317 m/s，相比試驗數據低約47%。綜合考慮，PR 方程更適用于分析高壓氫氣驅動的激波風洞流動問題，后續數值模擬采用PR 方程。

2.2 激波風洞流場結構和試驗氣流參數

運用前述準一維數值模擬方法對激波風洞的全場內流進行計算，通過對比采用立方型氣體狀態方程和理想氣體狀態方程的算例，考察高壓真實氣體效應導致的流場結構和流動參數的差異，考慮到試驗時間最大化的需求，本文全部工況中均在理想氣體狀態方程下通過調整驅動段和被驅動段氣體狀態，使激波風洞處于縫合接觸面運行狀態。同一工況、不同狀態方程的計算，激波管初始條件相同。

為避免引入多余擾動，計算過程考慮一個理想的瞬態、無損耗的破膜過程。是否考慮非理想破膜對于本文中關注的高壓真實氣體效應作用規律和機理沒有本質影響。

2.2.1 激波風洞全場流動結構

以表2 所示的典型工況為例考察激波風洞全場準一維流動的時空結構，該工況以150 MPa氫氣驅動110 kPa 氮氣。在使用理想氣體狀態方程并考慮管壁沿程損耗情況下，該工況中經過激波壓縮后的試驗氣流總參數為：總壓100 MPa、總溫6 000 K。

數值模擬獲得理想氣體狀態方程下激波風洞全場壓強、溫度x-t 云圖分別如圖7～8 所示。由圖7～8 可以看到：主膜片破裂產生一道右行入射激波、一道接觸面和一簇左行中心稀疏波。入射激波到達噴管入口后，喉道前方膜片破裂，氣流進入噴管并在喉道發生壅塞，從而反射一道激波向上游運動。喉道前方氣流在反射激波作用后成為高溫高壓的駐室氣流。

駐室氣流的參數及其穩定性直接決定了噴管出口試驗氣流的狀態及有效試驗時間。如圖7 所示，在縫合接觸面運行狀態下，對駐室氣流產生影響的波系主要為反射稀疏波頭和接觸面。定義入射激波抵達駐室時間為t1，反射稀疏波頭抵達駐室的時間為t2，接觸面抵達駐室的時間為t3，則駐室氣流參數相對穩定的有效試驗時間teff 可以定義為：

teff = min{t2 -t1，"t3 -t1} （20）

采用與表2 中相同的初始條件，使用PR 狀態方程進行考慮高壓真實氣體效應的模擬計算，得到的激波風洞全場的壓強、溫度x-t 云圖分別如圖9～10 所示。

對比2 種狀態方程下的整體流場結構。首先，兩者波系的時空結構特征相似，并且，入射激波以及接觸面的運動軌跡在激波反射前沒有明顯區別，入射激波均在7 ms 左右到達噴管入口并發生反射，說明入射激波強度相當；其次，激波反射之后，反射激波以及接觸面的后續運動狀態兩者存在可見差異，圖7中反射激波的運動速度明顯慢于圖9 中的，接觸面在與反射激波作用后，速度減小，緩慢向噴管移動，其中，圖7 中接觸面在約23 ms 進入噴管，而圖9 中接觸面在25 ms 內未能進入噴管；其三，稀疏波運動軌跡存在顯著差異，PR 方程下稀疏波和反射稀疏波速度顯著高于理想氣體，圖7 中稀疏波頭在約7.0 ms 時刻到達風洞左壁面，圖9 中稀疏波頭到達左壁面時間提前至約3.5 ms，反射稀疏波到達駐室的時間也大幅提前。

為了定量分析特征波系的運動情況，提取入射激波、反射激波、稀疏波頭和反射稀疏波頭運動速度的沿程分布情況，如圖11 所示。圖11 中空心點為IG 方程計算結果，實心點為PR 氣體方程計算結果。根據一維非定常氣動理論，稀疏波相對氣體以聲速傳播，左行稀疏波傳播速度為u?c，右行稀疏波傳播速度為u+c。由于驅動段氣體狀態均勻，因此左行中心稀疏波頭大致以恒定的聲速向風洞左端運動。其中，理想氣體方程下稀疏波頭運動速度為1 260 m/s，PR 氣體方程下稀疏波頭運動速度為2 370 m/s。中心稀疏波在管壁反射后向右傳播，同樣，PR 氣體方程下的反射稀疏波頭絕對運動速度大于理想氣體方程下的稀疏波速。從圖11 還可以看到，盡管被驅動段初始流場均勻，但受管壁摩擦以及熱傳導的影響，入射激波沿程強度有所衰減，激波運動速度逐漸降低，使用PR 氣體方程計算得到的入射激波初始運動由2 631 m/s 逐漸減低至2 544 m/s，使用理想氣體方程計算得到的入射激波初始運動由2 751 m/s 逐漸減低至2 635 m/s。PR 氣體方程計算得到的入射激波速度略低于理想氣體方程計算得到的入射激波初始速度。而對于反射激波，兩者對比趨勢相反，PR 氣體方程的反射激波速度大于理想氣體方程的反射激波速度。

關于氣體的熱化學非平衡效應，計算顯示，在喉道上游的激波管流動中，流場溫度尚不足以導致顯著的氮分子離解（盡管駐室區溫度可達6 000 K，但高壓抑制了離解的發生）；同時激波后振動弛豫的尺度遠小于管徑尺度，因此，整個激波管流動近似處于振動平衡和化學凍結的狀態。后續討論中將不再單獨討論振動溫度和離解組分。

2.2.2 駐室狀態

激波風洞內波系速度的變化會導致波系時空相干關系的改變，并共同導致駐室狀態參數和有效試驗時間的變化。圖12 為不考慮真實氣體效應（即，熱化學凍結（thermochemically frozen， TF），采用理想氣體狀態方程（IG EOS））、僅考慮高溫氣體效應（即，熱化學非平衡（thermochemical nonequilibrium，TN，采用理想氣體狀態方程（ IG EOS） ）和考慮高溫、高壓氣體效應（即，熱化學非平衡（ thermochemicalnonequilibrium，TN，采用PR 狀態方程（PR EOS））3 種情況下激波風洞駐室區（激波管右端，距離喉道0.15 m 處）的壓強、平動溫度和當地氣流組分摩爾分數隨時間變化情況。

如圖12 所示，主膜片破膜約7 ms 后，入射激波首先抵達右端，當地壓強、溫度躍升；在經歷短暫平臺后，反射激波抵達，壓強和溫度再次躍升；在經歷一段波動后，壓強和溫度信號趨于平穩，其中，壓強信號緩慢上升，溫度信號則呈緩慢降低趨勢，這一區域內參數無明顯波動，是激波風洞試驗主要利用的區域。

對比圖12 中3 種工況計算得到的駐室參數情況，激波管流動的高溫氣體效應主要體現為氮氣分子振動能激發和小部分氮分子的離解降低了駐室區的溫度。該效應對全場流動時空特性影響不大。為了分析高壓氣體效應的影響，后續不再詳述高溫氣體效應的影響，僅在同等考慮高溫氣體效應情況下，針對使用理想氣體狀態方程和PR 狀態方程的流動差異進行討論。

理想氣體狀態方程和PR 狀態方程下的駐室區壓強信號分別在20.6 和13.3 ms 時刻開始下降。對于溫度，理想氣體狀態方程下的駐室區溫度在16.7 ms 附近迅速下降；而PR 狀態方程下的溫度曲線無明顯趨勢改變，僅在13.3 ms 附近輕微向下折轉。由圖12 的組分變化可知，理想氣體下16.7 ms 后氫氣組分開始占據駐室區，表明此時接觸面抵達；而PR 方程下，25 ms 內駐室區均以氮氣為主，未見明顯的氫組分，表明接觸面尚未進入觀測區間。

結合波系結構和組分演變可知，圖12 中壓強的迅速下降主要對應于反射稀疏波的抵達，而溫度的迅速下降主要對應于接觸面（冷的驅動氣體）的抵達。由于試驗氣流需要同時維持壓強、溫度和組分的穩定，有效試驗時間由最先的抵達的稀疏波或接觸面決定。

提取各個波系到達駐室的時刻如表3 所示。受高壓真實氣體效應影響，入射激波抵達時間由6.8 ms 輕微延遲至7.1 ms，而反射稀疏波頭到達駐室的時間則由20.6 ms 大幅提前至13.5 ms，接觸面到達駐室時間由17.2 ms 延遲為25.5 ms。這使得限制有效試驗時間的因素由接觸面變為反射稀疏波，實際有效試驗時間由10.4 ms 縮短為6.4 ms，降幅達38.5%。

考察有效試驗時間內的駐室氣流參數情況。首先，由于PR 方程下入射激波強度較弱，因此其有效試驗時間內的壓強和溫度均略低于理想氣體計算結果。其次，由于入射激波沿程傳播受管壁摩擦和傳熱的影響強度逐步減弱，波后氣流非均勻，反射激波后的氣流參數不能維持平臺（試驗數據同樣顯示上述特征，如圖3 所示）。從圖12 可以看到，有效試驗時間內：在理想氣體方程下，駐室壓強由80 MPa逐步上升到120 MPa，平均駐室壓強約為100 MPa，駐室溫度由6 600 K 逐漸降低到6 100 K，平均駐室溫度約為6 350 K；而在PR 氣體方程下，駐室壓強由80 MPa 逐漸上升到100 MPa，平均駐室壓強約為90 MPa，駐室溫度由6 300 K 逐漸降低到5 800 K，平均駐室溫度約為6 050 K。

2.3 高壓真實氣體對激波管流動的影響機制

上述數值結果表明，考慮高壓真實氣體效應后風洞流場結構、試驗氣流參數和有效試驗時間均發生了明顯變化。為了揭示高壓真實氣體對激波管流動的影響機制，接下來將對采用理想氣體狀態方程和PR 狀態方程下激波管內稀疏波、激波以及后續接觸面的傳播過程進行對比分析。

2.3.1 稀疏波傳播速度

考慮高壓真實氣體效應后稀疏波以及反射稀疏波的傳播速度明顯加快。稀疏波以聲速傳播，因此稀疏波的傳播速度加快實際上是當地聲速在高壓真實氣體下發生增大所導致的。

根據式（19），以PR 方程為狀態方程時，氣體的聲速可寫作：

式中：根號下γ RT 為理想氣體狀態方程下的聲速。與狀態方程本身相對應，PR 方程下高壓真實氣體效應對聲速的影響可分解為分子體積和分子間作用力兩因素，分別由兩無量綱因子ηv 和ηF 表達。當分子體積與分子間作用力趨于0 時，ηv 和ηF 均趨于1，聲速回歸至理想氣體狀態方程下的聲速。

為了評估上述兩因素的影響程度，分別提取ηv 和ηF 隨壓強變化情況，如圖13 所示。可以看到，分子體積項的影響隨著壓強的增高而顯著增強，而分子間作用力的影響隨壓強增高變化不大，數值維持在1～0.95 之間。由此可知，分子體積對空間的占據效應對于聲速的影響占據主導地位。

2.3.2 激波管入射波系和驅動能力

除改變稀疏波傳播速度外，氣體狀態方程對激波和稀疏波的強度也存在不同程度的影響，這些影響可改變激波管和激波風洞的工作狀態。這里首先從氣體動力學基本方程出發分析真實氣體對激波管入射波系和驅動能力的影響機制。

激波管特征波系由向低壓區傳播的入射激波、向高壓區傳播的稀疏波和中間隔離驅動氣體與被驅動氣體的接觸面構成。三波系將激波管流動分為4 個平臺區，依次標記為區域①～④，入射激波在端壁發生反射將區域②氣流的動能轉化為氣體內能并疊加到區域②氣流既有的內能中，由此產生高溫和高壓的區域⑤。如圖14 所示為激波管流動波系和流場分區的示意圖。

由于高溫熱力學參數和PR 方程形式的復雜性，無法解析獲得高溫高壓真實氣體下的激波關系和稀疏波關系的具體形式。為此，只能由原始方程出發解算激波管的入射波系及其流場參數：（1） 入射激波關系由跨間斷的質量、動量、能量守恒方程外加PR 方程組成，4 個方程可建立波后壓強p、比容v、溫度T 和速度u 等4 個未知數和入射激波強度的關系。PR 方程和平衡熱力學參數主要通過能量方程中的焓h 發生作用。間斷關系方程全部為代數方程，因而以牛頓迭代數值求解（由已知的區域①參數出發，獲得區域②的參數）方程組：

式中：下標為1 的表示波前氣體參數，下標為2 的表示波后氣體參數。

（2） 稀疏波關系由小擾動波的連續性關系和等熵關系外加PR 方程組成，3 個方程可建立壓強、比容、溫度和速度（p、v、T、u）4 個未知數中任意兩者的關系。稀疏波關系中的相容關系為常微分方程，因而以R-K 方法積分求解（以區域④的條件為初始條件求解區域③的參數）如下聯立方程組：

式中：s、cV 分別為氣體的比熵、比定容熱容，cV，IG 為理想氣體狀態下氣體的比定壓熱容，c 和u 為當地聲速和當地速度。

區域②和區域③由接觸面隔開，因此區域②與區域③的速度（u）和壓強（p）應當分別相等。可在pu圖上尋找激波管問題的解點。

在激波管區域④和區域①參數已知的情況下，可分別繪制出入射激波后和稀疏波后氣流的p-u 圖，如圖15 所示。由圖15 可知，高壓真實氣體效應對稀疏波關系產生了較明顯的影響，驅動段氣體由相同壓強膨脹時，稀疏波后相同速度處基于PR 狀態方程的流場壓強更低，并且這種差異會隨著稀疏波的增強而擴大。與此同時，高壓真實氣體效應對入射激波間斷關系幾乎沒有影響，兩者的波后p-u 曲線幾乎重合（圖15 中灰色實線與紅色虛線）。這主要是因為：盡管激波后壓強顯著增大，但同時發生的溫度提升大幅壓制了氣體密度或分子數密度的增長，使得高壓氣體效應無法顯現。為了更清晰地體現這一效應，圖16 給出了PR 方程所定義的氮氣可壓縮性因子隨溫度和壓強的變化曲面（溫度300～4 000 K，壓強0.1～200 MPa），任意熱力學過程均對應于該曲面上的一條曲線。圖16 同時給出了初始溫度300 K、初始壓強0.1 MPa 情況下對應等溫、等熵和激波壓縮3 種路徑的參數變化曲線。可以看到：等溫壓縮路徑下溫度維持不變，可壓縮因子具有最顯著的上升趨勢；等熵壓縮路徑下溫度和可壓縮因子變化趨勢均居中；激波壓縮路徑下溫度上升最快，可壓縮因子則上升不到0.5%，表明該初始條件下的激波壓縮過程幾乎不呈現高壓真實氣體效應。

當激波管中激波和稀疏波后氣流在接觸面處匹配流場速度和壓強時（即圖15 中激波和稀疏波曲線交點），盡管激波關系幾乎不變，但由于稀疏波關系曲線的偏移，考慮高壓真實氣體效應獲得的激波強度仍要低于理想氣體情形。然而這里激波強度降低的幅度是有限的，理論計算顯示，在表1 條件下，激波管產生的初始激波馬赫數大致由7.9 降至7.6，這與準一維數值模擬結果是一致的。

為呈現廣泛條件下驅動段、被驅動段間壓強配比與入射激波間的關系，按照以上理論方法對2 種狀態方程下不同高、低壓條件驅動產生的入射激波強度進行計算。圖17 所示為以壓強200 和50 MPa 氫氣驅動壓強10～1 000 kPa 的氮氣的入射激波馬赫數，其中空心點為理想氣體狀態方程的計算結果，實心點為PR 狀態方程的計算結果。由圖17 可以看到：首先，考慮高壓效應的入射激波馬赫數要低于理想氣體情形，這一規律在全范圍內維持不變；其次，高壓真實氣體效應對入射激波馬赫數的影響幅度隨驅動段壓強和高、低壓段壓比的增大而增大，但整體變化幅度并不十分顯著，在所測試的最極端情形下（200 MPa 氫氣驅動10 kPa 氮氣，壓比20 000），理想氣體狀態方程下入射激波馬赫數約為11.37，PR 方程下激波馬赫數約為10.78，兩者差異僅為0.55。

2.3.3 激波管后續波系傳播

根據理論計算，考慮高壓真實氣體效應后，高壓驅動段的氣體聲速由理想氣體中1 317 m/s 提高至2 596 m/s，這使得以聲速傳播的入射稀疏波頭更快抵達左端壁。為進一步分析反射波系的傳播特征，將圖15 中理想氣體狀態方程和PR 狀態方程計算得到稀疏波后聲速和激波后聲速給出，如圖18 所示，圖18中2 條豎直線分別對應圖15 中高壓真實氣體和理想氣體的解。

由圖18 可知，真實氣體入射稀疏波后的驅動氣體聲速（圖18 中藍色實線）在高壓真實氣體解線之前始終高于理想氣體中相應部分的聲速（圖18 中黑色虛線）。因此，當真實氣體的左行稀疏波從端壁反射后，其傳播速度仍然顯著高于理想氣體。同時，真實氣體入射激波后聲速則略低于理想氣體；當右端反射激波與接觸面作用后，根據氣體動力學理論，由于真實氣體中接觸面左側（區域③）氣體有較高的聲速，反射激波進入左側氣體后有更快的傳播速度。

由于喉道處于壅塞狀態，反射激波作用后接觸面的運動速度很大程度上取決于駐室氣流的狀態。根據一維等熵流理論，壅塞流速大致正比于駐室區的聲速（或溫度0.5 次方），由于理想氣體中區域⑤的氣流溫度略高于真實氣體的（圖12），因此真實氣體接觸面的右行速度本來就略低于理想氣體，但兩者差異并不明顯。實際流動中真正使得真實氣體中接觸面延遲抵達的原因是反射稀疏波提前抵達并作用于駐室區，使得當地聲速急劇降低，從而極大降低了接觸面的后續運動速度。

2.4 不同條件下高壓效應對試驗時間的影響

激波風洞試驗往往需要在同一設備上獲取不同總溫和不同總壓的試驗氣流。通過改變激波風洞運行條件，獲得試驗氣流總溫2000～6 000 K、總壓20～200 MPa 范圍內縫合接觸面運行的不同試驗工況，以此考察不同條件下高壓真實氣體效應對試驗時間的影響。試驗氣流總溫的調整和縫合條件的實現通過調整驅動段組分（氫氣與氮氣的配比）和高、低壓段壓比實現，試驗氣流總壓的改變通過在維持高、低壓段壓比情況下調整壓強實現。各段初始氣體溫度恒為300 K。

首先，在維持試驗氣流名義總溫為6 000 K 的情況下，改變試驗氣流的總壓，分別提取理想氣體狀態方程和PR 狀態方程下數值模擬獲得的有效試驗時間進行比較。不同總壓條件下，入射激波、反射稀疏波頭和接觸面到達駐室的時刻（分別標記為t1、t2 和t3）如圖19 所示。

在理想氣體狀態方程下，由于激波管初始壓比和溫度維持不變，初始入射激波和稀疏波的速度均不變，而在湍流區激波沿程衰減幅度同樣變化不大，因此，各波系的時空關系基本維持不變，反射激波、反射稀疏波和接觸面的抵達時刻均無顯著改變。此時各工況均為接觸面率先抵達，有效試驗由t2?t1 決定。

在采用PR 氣體狀態方程之后，如前面的分析：隨著總壓的提高（驅動壓強增高），入射激波的強度基本未受影響（輕微減弱）；稀疏波傳播速度加快，反射稀疏波頭抵達時間提前；相應的，受提前抵達的反射稀疏波作用，駐室區氣流聲速降低，其經由喉道流出的速度降低，這使得接觸面抵達時間延遲。從圖19 可以看到，在試驗氣流總壓20 MPa 運行條件下，反射稀疏波提前幅度不大，接觸面仍在稀疏波前抵達端部，因此，接觸面并未出現顯著延遲，此時有效試驗時間仍由t2?t1 決定。而當總壓提高至50 MPa，反射稀疏波越過接觸面率先抵達端部，此時試驗氣流部分受到反射稀疏波作用，導致接觸面相對理想氣體情況發生滯后。這一效應隨著試驗氣流總壓的提升愈發顯著。總壓50 MPa 以上的工況，真實氣體有效試驗時間轉而由t3?t1 決定。在試驗氣流總壓100、150 和200 MPa 等3 種工況下，高壓真實氣體效應導致有效試驗時間相對理想氣體分別減少約40%、50% 和60%。

進一步對試驗氣流總溫2 000～6 000 K、總壓20～200 MPa 范圍的其余工況進行計算。為達到縫合接觸面運行條件，不同試驗氣流名義總溫對應激波管驅動段摻入氮氣摩爾分數以及高、低壓段的壓比如表4 所示。

將以上測試范圍內高壓真實氣體效應對有效試驗時間的影響情況總結為圖20，其中Δ 定義為2 種狀態方程下有效試驗時間的差值相對理想氣體方程有效試驗時間的比：

分析固定總壓條件下，有效試驗時間隨總溫的變化情況可知，隨著總溫升高， Δ 逐漸減小。結合2.2 節與2.3 節中的分析，溫度的影響主要體現在試驗氣流參數上，對于試驗時間影響較小。Δ 隨總溫升高逐漸減小的原因為高溫對高壓氣體效應的影響具有抑制作用。隨著溫度的升高，高壓氣體效應引起的稀疏波運動速度變化幅度減小，在限制有效試驗時間的因素為反射稀疏波到達后，理想氣體狀態方程和PR 狀態方程下有效試驗時間的差異進一步減小。

如圖20 所示，不同驅動狀態下高壓氣體效應的影響程度不同：（1） 在試驗氣流總溫較低時，試驗時間主要受到稀疏波限制，高壓真實氣體效應致使有效試驗時間直接呈現縮短趨勢。（2） 當試驗氣流總溫提高至3 000 K 以上，理想氣體中接觸面先于反射稀疏波抵達端部，此時如總壓較低，高壓效應不足以改變兩者的抵達順序，則有效試驗時間幾乎不變；而當總壓提升，高壓效應導致反射稀疏波先于接觸面抵達端部，則有效試驗時間隨試驗氣流總壓的提高而逐步縮短。

綜上所述，對于高壓氫氣驅動的激波風洞內的全場流動時空結構和駐室區氣流參數的分析獲得如下認知：（1） 高壓真實氣體效應對激波管流動的影響主要體現為稀疏波傳播速度加快侵蝕有效試驗時間長度；（2） 同等壓強下，溫度越高則高壓氣體效應越弱。因此，為了減小高壓真實氣體效應的影響，工程中易于實現的方法有：①延長驅動段的長度，使得稀疏波需要經過更長的距離才能反射到達駐室，延長有效試驗時間；②對激波管驅動段進行加熱處理（同氣體介質、同等壓強下，高溫氣體驅動能力更強，氣體密度更小，可削弱高壓氣體效應）。

3 結 論

推導了立方型氣體狀態方程下氣體的熱力學參數（內能、比熱容、比熱比等）和聲速，并將上述參數運用于耦合高溫熱化學非平衡過程的準一維數值模擬，對考慮高壓真實氣體效應的高焓激波風洞流動結構以及試驗氣流參數進行了準一維數值計算和分析，獲得了如下主要結論。

（1） 采用PR 氣體狀態方程可以很好地刻畫氣體的高壓真實氣體效應，耦合該方程的準一維數值模擬可以準確復現高焓激波風洞的整體流動。

（2） 高溫真實氣體效應和高壓真實氣體效應在物理上近乎處于相互排斥的地位，兩者分別在激波管流動的熱端（被驅動氣體）和冷端（驅動氣體）起主導性作用。高壓真實氣體效應對激波管的驅動能力以及入射激波和反射激波后的流動狀態影響甚微，僅在驅動壓強較高時輕微削弱所產生激波的強度；但在溫度不高的驅動氣體中，高壓真實氣體效應使氣體聲速增大，從而導致稀疏波和反射稀疏波傳播速度加快，這一效應可以極大地改變激波管中各波系的相干時空關系。

（3） 高壓真實氣體效應對激波風洞駐室試驗氣流的影響主要呈現為它對有效試驗時間長度的侵蝕。這種作用主要通過高壓效應加快稀疏波的傳播、從而致使稀疏波提前抵達駐室區并顯著降低當地氣流總壓達成。如高壓效應不能使反射稀疏波先于接觸面抵達喉道，則高壓效應不對有效試驗時間形成實質侵蝕。

（4） 激波風洞在試驗氣流總溫較低而總壓較高的工況下，有效試驗時間一般以反射稀疏波頭抵達時刻為終點，高壓真實氣體效應更容易在此類工況下對有效試驗時間產生大幅影響。對于采用冷高壓氣體驅動的高焓激波風洞，可在理想設計的基礎上適當延長驅動段長度以抵消反射稀疏波提前抵達對試驗時間的影響。

參考文獻：

[1]ANDERSON J D JR. Hypersonic and high-temperature gas dynamics [M]. 2nd ed. Reston， Virginia， USA： AIAA， 2006： 449–461. DOI： 10.2514/4.861956.

[2] VAN DER WAALS J D. On the continuity of the gaseous and liquid state [D]. Leiden， Holland： University of Leiden， 1873.

[3]JACOBS P A. Quasi-one-dimensional modeling of a free-piston shock tunnel [J]. AIAA Journal， 1994， 32（1）： 137–145. DOI：10.2514/3.11961.

[4]BOYCE R R， TAKAHASHI M， STALKER R J. Mass spectrometric measurements of driver gas arrival in the T4 free-piston"shock-tunnel [J]. Shock Waves， 2005， 14（5/6）： 371–378. DOI： 10.1007/s00193-005-0276-3.

[5]HANNEMANN K. High enthalpy flows in the HEG shock tunnel： experiment and numerical rebuilding [C]//41st Aerospace"Sciences Meeting and Exhibit. Reno： AIAA， 2003： 978. DOI： 10.2514/6.2003-978.

[6]JOHNSTON I， WEILAND M， SCHRAMM J， et al. Aerothermodynamics of the ARD-Postflight numerics and shock-tunnel"experiments [C]//40th AIAA Aerospace Sciences Meeting amp; Exhibit. Reno： AIAA， 2002： 407. DOI： 10.2514/6.2002-407.

[7]HOLDEN M， WADHAMS T， MACLEAN M， et al. Experimental studies in LENS I and X to evaluate real gas effects on"hypevelocity vehicle performance [C]//45th AIAA Aerospace Sciences Meeting and Exhibit. Reno： AIAA， 2007： 204. DOI：10.2514/6.2007-204.

[8]盧洪波， 陳星， 曾憲政， 等. FD-21 風洞Ma=10 高超聲速推進試驗技術探索 [J]. 氣體物理， 2022， 7（2）： 1–12. DOI： 10.19527/j.cnki.2096-1642.0926.

LU H B， CHEN X， ZENG X Z， et al. Exploration of experimental techniques on Ma=10 scramjets in FD-21 high enthalpy"shock tunnel [J]. Physics of Gases， 2022， 7（2）： 1–12. DOI： 10.19527/j.cnki.2096-1642.0926.

[9]LLOBET J R， YAMADA G， TONIATO P. Quasi-0D-model mapping of reflected shock tunnel operating conditions [J].AIAA Journal， 2020， 58（4）： 1668–1680. DOI： 10.2514/1.J058660.

[10]LYNCH K P， GRASSER T， FARIAS P， et al. Design and characterization of the Sandia free-piston reflected shock tunnel [C]//AIAA SCITECH 2022 Forum. San Diego： AIAA， 2022： 0968. DOI： 10.2514/6.2022-0968.

[11]GORDON S， MCBRIDE B J. Computer program for calculation of complex chemical equilibrium compositions and"applications. Part 1： analysis： NASA RP-1311 [R]. Cleveland， Ohio， USA： NASA， 1996.

[12]LUO K， WANG Q， LI J W， et al. Numerical modeling of a high-enthalpy shock tunnel driven by gaseous detonation [J].Aerospace Science and Technology， 2020， 104： 105958. DOI： 10.1016/j.ast.2020.105958.

[13]MACLEAN M， WADHAMS T， HOLDEN M， et al. Investigation of blunt bodies with CO2 test gas including catalytic effects [C]//38th AIAA Thermophysics Conference. Toronto： AIAA， 2005： 4693. DOI： 10.2514/6.2005-4693.

[14]MACLEAN M， DUFRENE A， WADHAMS T， et al. Numerical and experimental characterization of high enthalpy flow in an"expansion tunnel facility [C]//48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace"Exposition. Orlando： AIAA， 2010： 1562. DOI： 10.2514/6.2010-1562.

[15]CANDLER G. Hypersonic nozzle analysis using an excluded volume equation of state [C]//38th AIAA Thermophysics"Conference. Toronto： AIAA， 2005： 5202. DOI： 10.2514/6.2005-5202.

[16]SIRIGNANO W A. Compressible flow at high pressure with linear equation of state [J]. Journal of Fluid Mechanics， 2018，843： 244–292. DOI： 10.1017/jfm.2018.166.

[17]王剛， 馬曉偉， 江濤， 等. 前向空腔的弓形激波特性研究 [J]. 宇航學報， 2016， 37（8）： 901–907. DOI： 10.3873/j.issn.1000-1328.2016.08.002.

WANG G， MA X W， JIANG T， et al. Characteristics of bow shock for a forward-facing cavity [J]. Journal of Astronautics，2016， 37（8）： 901–907. DOI： 10.3873/j.issn.1000-1328.2016.08.002.

[18]GROTH C P T， GOTTLIEB J J. Numerical study of two-stage light-gas hypervelocity projectile launchers： UTIAS Report"372 [R]. Toronto： University of Toronto Intsitute for Aerospace Studies， 1988.

[19]WILKE C R. A viscosity equation for gas mixtures [J]. The Journal of Chemical Physics， 1950， 18（4）： 517–519. DOI： 10.1063/1.1747673.

[20]PARK C. Assessment of two-temperature kinetic model for ionizing air [J]. Journal of Thermophysics and Heat Transfer，1989， 3（3）： 233–244. DOI： 10.2514/3.28771.

[21]MILLIKAN R C， WHITE D R. Systematics of vibrational relaxation [J]. The Journal of Chemical Physics， 1963， 39（12）：3209–3213. DOI： 10.1063/1.1734182.

[22]GUPTA R N， YOS J M， THOMPSON R A. A review of reaction rates and thermodynamic and transport properties for the 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K： NASA TM-101528 [R]. Hampton：NASA， 1990.

[23]DUNN M G， KANG S. Theoretical and experimental studies of reentry plasmas： NASA CR-2232 [R]. Buffalo， New York，USA： NASA， 1973.

[24]TORO E F， SPRUCE M， SPEARES W. Restoration of the contact surface in the HLL-Riemann solver [J]. Shock Waves，1994， 4（1）： 25–34. DOI： 10.1007/bf01414629.

[25]侯自豪， 朱雨建， 薄靖龍， 等. 真空管道列車準一維氣動特性 [J]. 機械工程學報， 2022， 58（6）： 119–129. DOI：10.3901/JME.2022.06.119.

HOU Z H， ZHU Y J， BO J L， et al. Quasi-one-dimensional aerodynamic characteristics of tube train [J]. Journal of Mechanical"Engineering， 2022， 58（6）： 119–129. DOI： 10.3901/JME.2022.06.119.

[26]HOU Z H， ZHU Y J， BO J L， et al. A quasi-one-dimensional study on global characteristics of tube train flows [J]. Physics of Fluids， 2022， 34（2）： 026104. DOI： 10.1063/5.0080544.

[27]YANG J T， ZHU Y J， YANG J M. Self-ignition induced by cylindrically imploding shock adapting to a convergent channel [J].Physics of Fluids， 2017， 29（3）： 031702. DOI： 10.1063/1.4979135.

[28]中國空氣動力研究與發展中心超高速空氣動力研究所， 中國科學技術大學. 超高壓激波類設備熱化學非平衡準一維計算軟件. V1.0： 2022SR0066696 [P]. 2022-01-11.

[29]PENG D Y， ROBINSON D B. A new two-constant equation of state [J]. Industrial and Engineering Chemistry Fundamentals，1976， 15（1）： 59–64. DOI： 10.1021/i160057a011.

[30]LOPEZ-ECHEVERRY J S， REIF-ACHERMAN S， ARAUJO-LOPEZ E. Peng-Robinson equation of state： 40 years through"cubics [J]. Fluid Phase Equilibria， 2017， 447： 39–71. DOI： 10.1016/j.fluid.2017.05.007.

[31]SOAVE G. Equilibrium constants from a modified Redlich-Kwong equation of state [J]. Chemical Engineering Science， 1972，27（6）： 1197–1203. DOI： 10.1016/0009-2509（72）80096-4.

[32]ZUDKEVITCH D， JOFFE J. Correlation and prediction of vapor-liquid equilibria with the Redlich-Kwong equation of"state [J]. AIChE Journal， 1970， 16（1）： 112–119. DOI： 10.1002/aic.690160122.

[33]POLING B E， PRAUSNITZ J M， O’CONNELL J P. Properties of gases and liquids [M]. 5th ed. New York： McGraw-Hill"Education， 2001： 151-155.

[34]KAY W B. Gases and vapors at high temperature and pressure： density of hydrocarbon [J]. Industrial and Engineering"Chemistry， 1936， 28（9）： 1014–1019. DOI： 10.1021/ie50321a008.

[35]LINSTROM P J， MALLARD W G. The NIST Chemistry WebBook： a chemical data resource on the internet [J]. Journal of"Chemical and Engineering Data， 2001， 46（5）： 1059–1063. DOI： 10.1021/je000236i.

（責任編輯 張凌云）

基金項目： 國家自然科學基金（12172354，U21B6003）