基于立方型氣體狀態(tài)方程的激波風(fēng)洞準(zhǔn)一維流動(dòng)數(shù)值研究

摘要： 采用基于立方型氣體狀態(tài)方程的準(zhǔn)一維流動(dòng)數(shù)值模擬方法研究了反射式高焓激波風(fēng)洞的真實(shí)氣體流動(dòng)，重點(diǎn)關(guān)注了高壓真實(shí)氣體效應(yīng)對(duì)風(fēng)洞全場(chǎng)流動(dòng)時(shí)空結(jié)構(gòu)和駐室區(qū)氣流參數(shù)的影響，并以理論分析揭示了高壓真實(shí)氣體效應(yīng)對(duì)激波管內(nèi)流動(dòng)的作用機(jī)理。研究表明：對(duì)于以冷高壓氣體驅(qū)動(dòng)的激波風(fēng)洞，使用考慮分子體積和分子間作用力的真實(shí)氣體狀態(tài)方程能夠更準(zhǔn)確地描述氣體的狀態(tài)和風(fēng)洞內(nèi)的流動(dòng)狀況。高壓真實(shí)氣體效應(yīng)主要在冷驅(qū)動(dòng)氣體中發(fā)生作用，其作用效果主要是使當(dāng)?shù)芈曀僭龃螅瑥亩沟萌肷湎∈璨ê头瓷湎∈璨ǖ膫鞑ニ俣燃涌欤涣硪环矫妫邏簹怏w效應(yīng)在高溫氣體效應(yīng)較顯著的被驅(qū)動(dòng)氣體中作用微弱，且對(duì)激波管產(chǎn)生激波的強(qiáng)度和激波后的流動(dòng)狀態(tài)影響甚微。稀疏波的加快傳播改變了激波管波系的相干時(shí)空關(guān)系。提前抵達(dá)的稀疏波可在一定情況下侵蝕激波風(fēng)洞的有效試驗(yàn)時(shí)間。對(duì)于所測(cè)試的激波風(fēng)洞構(gòu)型，在150 MPa 氫氣驅(qū)動(dòng)110 kPa 氮?dú)獾墓r下，高壓效應(yīng)導(dǎo)致的有效試驗(yàn)時(shí)間縮短約38%。適當(dāng)加長(zhǎng)驅(qū)動(dòng)段長(zhǎng)度和采用高溫氣體驅(qū)動(dòng)均可有效減弱高壓真實(shí)氣體效應(yīng)的影響。

關(guān)鍵詞： 激波風(fēng)洞；高壓真實(shí)氣體效應(yīng)；立方型氣體狀態(tài)方程；準(zhǔn)一維數(shù)值模擬；有效試驗(yàn)時(shí)間

中圖分類號(hào)： O354.7 國(guó)標(biāo)學(xué)科代碼： 13025 文獻(xiàn)標(biāo)志碼： A

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

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

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

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

1 風(fēng)洞構(gòu)型與數(shù)值模擬方法

圖1 為FD-14A 激波風(fēng)洞[17] 內(nèi)流道的基本構(gòu)型。該風(fēng)洞激波管部分管道的內(nèi)徑均為0.15 m；當(dāng)前研究設(shè)置驅(qū)動(dòng)段長(zhǎng)9 m，被驅(qū)動(dòng)段長(zhǎng)18 m。為考慮反射端的質(zhì)量通量，設(shè)置噴管喉道的直徑為23.9 mm。

確定喉道和收縮段后，噴管其余部分不會(huì)對(duì)實(shí)際激波管流動(dòng)產(chǎn)生顯著影響，為了保留激波風(fēng)洞的全部結(jié)構(gòu)特征，計(jì)算過(guò)程中保留了整個(gè)噴管。

該風(fēng)洞驅(qū)動(dòng)段氣體為高壓氫氣（H2），被驅(qū)動(dòng)段氣體為空氣或氮?dú)猓∟2）。本研究中，為重點(diǎn)考察高壓真實(shí)氣體效應(yīng)，被驅(qū)動(dòng)段統(tǒng)一采用氮?dú)狻?/p>

1.1 準(zhǔn)一維流動(dòng)控制方程

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

式中：U 為守恒項(xiàng)；F 為對(duì)流項(xiàng)；P 為變截面管道壁面的流向作用力；A 為流道面積函數(shù)，描述管道截面空間變化；ρ、u、p、e、ev、Yi 分別為氣體的密度、速度、壓強(qiáng)、質(zhì)量總能量、振動(dòng)能和組分 的質(zhì)量分?jǐn)?shù)。非平衡氣體的質(zhì)量總能量有：

式中：Ru 為普適氣體常數(shù)，wi 為組分i 的分子摩爾質(zhì)量，a1，i、a2，i、a3，i、a4，i 和a5，i 為該組分的熱力學(xué)系數(shù)。Δep 為高壓真實(shí)氣體效應(yīng)的能量修正項(xiàng)，與真實(shí)氣體狀態(tài)方程相關(guān)，后續(xù)章節(jié)將給出推導(dǎo)形式。

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

式中：λi 為由組分?jǐn)U散導(dǎo)致的第i 組分含量變化，τ 為由黏性引起的流體內(nèi)部應(yīng)力，qtr 和qv 分別表示平動(dòng)能以及振動(dòng)能的熱傳導(dǎo)，M 表示壁面黏性導(dǎo)致的動(dòng)量傳遞，Qtr 和Qv 分別表示壁面?zhèn)鳠釋?dǎo)致的平動(dòng)能和振動(dòng)能傳遞， 表示由化學(xué)反應(yīng)導(dǎo)致的第i 組分生成，sv 表示由振動(dòng)弛豫導(dǎo)致的能量轉(zhuǎn)化。對(duì)于計(jì)算過(guò)程中需使用的氣體輸運(yùn)系數(shù)（主要是黏性系數(shù)和熱傳導(dǎo)系數(shù)），單一組分的輸運(yùn)系數(shù)以Sutherland 公式計(jì)算，混合氣體輸運(yùn)系數(shù)由單一組分系數(shù)按Wilke[19] 的方法混合得到。

壁面動(dòng)量損失源項(xiàng)M 和壁面能量損失源項(xiàng)Q 的表達(dá)式分別為：

式中：D 為管道水力學(xué)直徑；f 為Darcy-Weisbach 摩擦因數(shù)[18]；cp 為流體的比定壓熱容；Pr 為普朗特?cái)?shù)；Tw 為當(dāng)?shù)乇跍兀籘aw 對(duì)平動(dòng)能為絕熱壁溫，取值為經(jīng)過(guò)可壓縮修正之后的氣流平均溫度，對(duì)振動(dòng)能為當(dāng)?shù)仄骄駝?dòng)溫度。

本文中振動(dòng)弛豫采用Park[20] 的雙溫度模型，特征弛豫時(shí)間由Millikan 等[21] 的半經(jīng)驗(yàn)表達(dá)式計(jì)算得到，并參考Park[20] 的方法對(duì)極高溫度下的弛豫時(shí)間進(jìn)行修正。空氣在高溫下的化學(xué)反應(yīng)采用Gupta 等[22]的基元反應(yīng)動(dòng)力學(xué)機(jī)理（本文中僅涉及氫氣和氮?dú)猓蚨鴥H取其中氮?dú)怆x解反應(yīng)）。具體模型和數(shù)據(jù)可參考文獻(xiàn)[22-23]。

為了在統(tǒng)一模型架構(gòu)下求解激波風(fēng)洞的高溫高壓真實(shí)氣體流動(dòng)，在考慮高溫氣體效應(yīng)的基礎(chǔ)上，采用以下立方型氣體狀態(tài)方程：

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

1.2 數(shù)值計(jì)算方法

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

1.3 數(shù)值方法驗(yàn)證

本研究基于自主開發(fā)的準(zhǔn)一維熱化學(xué)非平衡數(shù)值模擬程序，所采用的基本模型和方法已在1.1 節(jié)和1.2 節(jié)中描述。該程序的基本算法和輸運(yùn)模型已在多個(gè)研究中獲得充分驗(yàn)證[25-27]，并在實(shí)際脈沖風(fēng)洞結(jié)構(gòu)和工況設(shè)計(jì)中獲得應(yīng)用[28]。

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

圖2 顯示，準(zhǔn)一維數(shù)值模擬可很好地刻畫廣泛條件下運(yùn)動(dòng)激波由于壁面摩擦和傳熱而發(fā)生衰減的現(xiàn)象。

接下來(lái)，對(duì)引入立方型狀態(tài)方程后的程序進(jìn)行補(bǔ)充驗(yàn)證。選取表1 所示的典型工況，計(jì)算激波風(fēng)洞駐室壓強(qiáng)隨時(shí)間演化情況，并與試驗(yàn)測(cè)量壓強(qiáng)數(shù)據(jù)進(jìn)行對(duì)比，如圖3 所示。

由圖3 可知，數(shù)值模擬計(jì)算得到激波到達(dá)時(shí)刻以及駐室區(qū)壓強(qiáng)變化情況與試驗(yàn)數(shù)據(jù)基本吻合。考慮立方型氣體狀態(tài)方程后的數(shù)值模擬更準(zhǔn)確地捕捉到了稀疏波系的運(yùn)動(dòng)情況，計(jì)算結(jié)果與試驗(yàn)結(jié)果的吻合度顯著提升。以上結(jié)果表明，本文中采用的準(zhǔn)一維數(shù)值模擬方法和程序可以準(zhǔn)確預(yù)測(cè)激波風(fēng)洞的整體運(yùn)行狀態(tài)。

2 結(jié)果與討論

2.1 立方型氣體狀態(tài)方程和熱力學(xué)參數(shù)

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

式中：參數(shù)b 表征分子體積，參數(shù)a 度量分子間吸引力；考慮到分子體積會(huì)對(duì)分子間作用力產(chǎn)生影響，因此引力項(xiàng)中還包含關(guān)于體積（比容v）的函數(shù)g（v），VDW 方程中取g（v）=v2。參數(shù)a、b 與氣體臨界參數(shù)（臨界溫度Tcr、臨界壓強(qiáng)pcr）相關(guān)，其中，a0 和b0 為狀態(tài)方程設(shè)定的常系數(shù)：

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

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

式中：fω 為分子偏心因子ω 的函數(shù)。

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

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

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

將2 種混合規(guī)則計(jì)算得到的70% 氫氣和30% 氮?dú)饣旌蠚怏w的可壓縮性因子（z=pv/RT）隨氣體壓強(qiáng)p 的變化與NIST 試驗(yàn)數(shù)據(jù)[35] 進(jìn)行對(duì)比，如圖4 所示。可以看到，使用經(jīng)典混合規(guī)則計(jì)算方法獲得的可壓縮性因子大幅偏離試驗(yàn)值，而使用準(zhǔn)臨界參數(shù)混合規(guī)則可以較好地描述混合氣體狀態(tài)。因此，后續(xù)計(jì)算將使用式（14） 進(jìn)行平均處理。

氣體的熱力學(xué)參數(shù)（內(nèi)能、焓、熵和比熱等）的具體形式與氣體狀態(tài)方程密切相關(guān)，氣體狀態(tài)方程變化后，需要對(duì)上述參數(shù)進(jìn)行重新推導(dǎo)。針對(duì)立方型方程一般形式（6），由熱力學(xué)基本方程推導(dǎo)獲得真實(shí)氣體狀態(tài)方程與理想氣體狀態(tài)方程內(nèi)能之差（式（2） 中Δep 項(xiàng)），表達(dá)如下：

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

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

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

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

2.2 激波風(fēng)洞流場(chǎng)結(jié)構(gòu)和試驗(yàn)氣流參數(shù)

運(yùn)用前述準(zhǔn)一維數(shù)值模擬方法對(duì)激波風(fēng)洞的全場(chǎng)內(nèi)流進(jìn)行計(jì)算，通過(guò)對(duì)比采用立方型氣體狀態(tài)方程和理想氣體狀態(tài)方程的算例，考察高壓真實(shí)氣體效應(yīng)導(dǎo)致的流場(chǎng)結(jié)構(gòu)和流動(dòng)參數(shù)的差異，考慮到試驗(yàn)時(shí)間最大化的需求，本文全部工況中均在理想氣體狀態(tài)方程下通過(guò)調(diào)整驅(qū)動(dòng)段和被驅(qū)動(dòng)段氣體狀態(tài)，使激波風(fēng)洞處于縫合接觸面運(yùn)行狀態(tài)。同一工況、不同狀態(tài)方程的計(jì)算，激波管初始條件相同。

為避免引入多余擾動(dòng)，計(jì)算過(guò)程考慮一個(gè)理想的瞬態(tài)、無(wú)損耗的破膜過(guò)程。是否考慮非理想破膜對(duì)于本文中關(guān)注的高壓真實(shí)氣體效應(yīng)作用規(guī)律和機(jī)理沒(méi)有本質(zhì)影響。

2.2.1 激波風(fēng)洞全場(chǎng)流動(dòng)結(jié)構(gòu)

以表2 所示的典型工況為例考察激波風(fēng)洞全場(chǎng)準(zhǔn)一維流動(dòng)的時(shí)空結(jié)構(gòu)，該工況以150 MPa氫氣驅(qū)動(dòng)110 kPa 氮?dú)狻Ｔ谑褂美硐霘怏w狀態(tài)方程并考慮管壁沿程損耗情況下，該工況中經(jīng)過(guò)激波壓縮后的試驗(yàn)氣流總參數(shù)為：總壓100 MPa、總溫6 000 K。

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

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

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

采用與表2 中相同的初始條件，使用PR 狀態(tài)方程進(jìn)行考慮高壓真實(shí)氣體效應(yīng)的模擬計(jì)算，得到的激波風(fēng)洞全場(chǎng)的壓強(qiáng)、溫度x-t 云圖分別如圖9～10 所示。

對(duì)比2 種狀態(tài)方程下的整體流場(chǎng)結(jié)構(gòu)。首先，兩者波系的時(shí)空結(jié)構(gòu)特征相似，并且，入射激波以及接觸面的運(yùn)動(dòng)軌跡在激波反射前沒(méi)有明顯區(qū)別，入射激波均在7 ms 左右到達(dá)噴管入口并發(fā)生反射，說(shuō)明入射激波強(qiáng)度相當(dāng)；其次，激波反射之后，反射激波以及接觸面的后續(xù)運(yùn)動(dòng)狀態(tài)兩者存在可見差異，圖7中反射激波的運(yùn)動(dòng)速度明顯慢于圖9 中的，接觸面在與反射激波作用后，速度減小，緩慢向噴管移動(dòng)，其中，圖7 中接觸面在約23 ms 進(jìn)入噴管，而圖9 中接觸面在25 ms 內(nèi)未能進(jìn)入噴管；其三，稀疏波運(yùn)動(dòng)軌跡存在顯著差異，PR 方程下稀疏波和反射稀疏波速度顯著高于理想氣體，圖7 中稀疏波頭在約7.0 ms 時(shí)刻到達(dá)風(fēng)洞左壁面，圖9 中稀疏波頭到達(dá)左壁面時(shí)間提前至約3.5 ms，反射稀疏波到達(dá)駐室的時(shí)間也大幅提前。

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

關(guān)于氣體的熱化學(xué)非平衡效應(yīng)，計(jì)算顯示，在喉道上游的激波管流動(dòng)中，流場(chǎng)溫度尚不足以導(dǎo)致顯著的氮分子離解（盡管駐室區(qū)溫度可達(dá)6 000 K，但高壓抑制了離解的發(fā)生）；同時(shí)激波后振動(dòng)弛豫的尺度遠(yuǎn)小于管徑尺度，因此，整個(gè)激波管流動(dòng)近似處于振動(dòng)平衡和化學(xué)凍結(jié)的狀態(tài)。后續(xù)討論中將不再單獨(dú)討論振動(dòng)溫度和離解組分。

2.2.2 駐室狀態(tài)

激波風(fēng)洞內(nèi)波系速度的變化會(huì)導(dǎo)致波系時(shí)空相干關(guān)系的改變，并共同導(dǎo)致駐室狀態(tài)參數(shù)和有效試驗(yàn)時(shí)間的變化。圖12 為不考慮真實(shí)氣體效應(yīng)（即，熱化學(xué)凍結(jié)（thermochemically frozen， TF），采用理想氣體狀態(tài)方程（IG EOS））、僅考慮高溫氣體效應(yīng)（即，熱化學(xué)非平衡（thermochemical nonequilibrium，TN，采用理想氣體狀態(tài)方程（ IG EOS） ）和考慮高溫、高壓氣體效應(yīng)（即，熱化學(xué)非平衡（ thermochemicalnonequilibrium，TN，采用PR 狀態(tài)方程（PR EOS））3 種情況下激波風(fēng)洞駐室區(qū)（激波管右端，距離喉道0.15 m 處）的壓強(qiáng)、平動(dòng)溫度和當(dāng)?shù)貧饬鹘M分摩爾分?jǐn)?shù)隨時(shí)間變化情況。

如圖12 所示，主膜片破膜約7 ms 后，入射激波首先抵達(dá)右端，當(dāng)?shù)貕簭?qiáng)、溫度躍升；在經(jīng)歷短暫平臺(tái)后，反射激波抵達(dá)，壓強(qiáng)和溫度再次躍升；在經(jīng)歷一段波動(dòng)后，壓強(qiáng)和溫度信號(hào)趨于平穩(wěn)，其中，壓強(qiáng)信號(hào)緩慢上升，溫度信號(hào)則呈緩慢降低趨勢(shì)，這一區(qū)域內(nèi)參數(shù)無(wú)明顯波動(dòng)，是激波風(fēng)洞試驗(yàn)主要利用的區(qū)域。

對(duì)比圖12 中3 種工況計(jì)算得到的駐室參數(shù)情況，激波管流動(dòng)的高溫氣體效應(yīng)主要體現(xiàn)為氮?dú)夥肿诱駝?dòng)能激發(fā)和小部分氮分子的離解降低了駐室區(qū)的溫度。該效應(yīng)對(duì)全場(chǎng)流動(dòng)時(shí)空特性影響不大。為了分析高壓氣體效應(yīng)的影響，后續(xù)不再詳述高溫氣體效應(yīng)的影響，僅在同等考慮高溫氣體效應(yīng)情況下，針對(duì)使用理想氣體狀態(tài)方程和PR 狀態(tài)方程的流動(dòng)差異進(jìn)行討論。

理想氣體狀態(tài)方程和PR 狀態(tài)方程下的駐室區(qū)壓強(qiáng)信號(hào)分別在20.6 和13.3 ms 時(shí)刻開始下降。對(duì)于溫度，理想氣體狀態(tài)方程下的駐室區(qū)溫度在16.7 ms 附近迅速下降；而PR 狀態(tài)方程下的溫度曲線無(wú)明顯趨勢(shì)改變，僅在13.3 ms 附近輕微向下折轉(zhuǎn)。由圖12 的組分變化可知，理想氣體下16.7 ms 后氫氣組分開始占據(jù)駐室區(qū)，表明此時(shí)接觸面抵達(dá)；而PR 方程下，25 ms 內(nèi)駐室區(qū)均以氮?dú)鉃橹鳎匆娒黠@的氫組分，表明接觸面尚未進(jìn)入觀測(cè)區(qū)間。

結(jié)合波系結(jié)構(gòu)和組分演變可知，圖12 中壓強(qiáng)的迅速下降主要對(duì)應(yīng)于反射稀疏波的抵達(dá)，而溫度的迅速下降主要對(duì)應(yīng)于接觸面（冷的驅(qū)動(dòng)氣體）的抵達(dá)。由于試驗(yàn)氣流需要同時(shí)維持壓強(qiáng)、溫度和組分的穩(wěn)定，有效試驗(yàn)時(shí)間由最先的抵達(dá)的稀疏波或接觸面決定。

提取各個(gè)波系到達(dá)駐室的時(shí)刻如表3 所示。受高壓真實(shí)氣體效應(yīng)影響，入射激波抵達(dá)時(shí)間由6.8 ms 輕微延遲至7.1 ms，而反射稀疏波頭到達(dá)駐室的時(shí)間則由20.6 ms 大幅提前至13.5 ms，接觸面到達(dá)駐室時(shí)間由17.2 ms 延遲為25.5 ms。這使得限制有效試驗(yàn)時(shí)間的因素由接觸面變?yōu)榉瓷湎∈璨ǎ瑢?shí)際有效試驗(yàn)時(shí)間由10.4 ms 縮短為6.4 ms，降幅達(dá)38.5%。

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

2.3 高壓真實(shí)氣體對(duì)激波管流動(dòng)的影響機(jī)制

上述數(shù)值結(jié)果表明，考慮高壓真實(shí)氣體效應(yīng)后風(fēng)洞流場(chǎng)結(jié)構(gòu)、試驗(yàn)氣流參數(shù)和有效試驗(yàn)時(shí)間均發(fā)生了明顯變化。為了揭示高壓真實(shí)氣體對(duì)激波管流動(dòng)的影響機(jī)制，接下來(lái)將對(duì)采用理想氣體狀態(tài)方程和PR 狀態(tài)方程下激波管內(nèi)稀疏波、激波以及后續(xù)接觸面的傳播過(guò)程進(jìn)行對(duì)比分析。

2.3.1 稀疏波傳播速度

考慮高壓真實(shí)氣體效應(yīng)后稀疏波以及反射稀疏波的傳播速度明顯加快。稀疏波以聲速傳播，因此稀疏波的傳播速度加快實(shí)際上是當(dāng)?shù)芈曀僭诟邏赫鎸?shí)氣體下發(fā)生增大所導(dǎo)致的。

根據(jù)式（19），以PR 方程為狀態(tài)方程時(shí)，氣體的聲速可寫作：

式中：根號(hào)下γ RT 為理想氣體狀態(tài)方程下的聲速。與狀態(tài)方程本身相對(duì)應(yīng)，PR 方程下高壓真實(shí)氣體效應(yīng)對(duì)聲速的影響可分解為分子體積和分子間作用力兩因素，分別由兩無(wú)量綱因子ηv 和ηF 表達(dá)。當(dāng)分子體積與分子間作用力趨于0 時(shí)，ηv 和ηF 均趨于1，聲速回歸至理想氣體狀態(tài)方程下的聲速。

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

2.3.2 激波管入射波系和驅(qū)動(dòng)能力

除改變稀疏波傳播速度外，氣體狀態(tài)方程對(duì)激波和稀疏波的強(qiáng)度也存在不同程度的影響，這些影響可改變激波管和激波風(fēng)洞的工作狀態(tài)。這里首先從氣體動(dòng)力學(xué)基本方程出發(fā)分析真實(shí)氣體對(duì)激波管入射波系和驅(qū)動(dòng)能力的影響機(jī)制。

激波管特征波系由向低壓區(qū)傳播的入射激波、向高壓區(qū)傳播的稀疏波和中間隔離驅(qū)動(dòng)氣體與被驅(qū)動(dòng)氣體的接觸面構(gòu)成。三波系將激波管流動(dòng)分為4 個(gè)平臺(tái)區(qū)，依次標(biāo)記為區(qū)域①～④，入射激波在端壁發(fā)生反射將區(qū)域②氣流的動(dòng)能轉(zhuǎn)化為氣體內(nèi)能并疊加到區(qū)域②氣流既有的內(nèi)能中，由此產(chǎn)生高溫和高壓的區(qū)域⑤。如圖14 所示為激波管流動(dòng)波系和流場(chǎng)分區(qū)的示意圖。

由于高溫?zé)崃W(xué)參數(shù)和PR 方程形式的復(fù)雜性，無(wú)法解析獲得高溫高壓真實(shí)氣體下的激波關(guān)系和稀疏波關(guān)系的具體形式。為此，只能由原始方程出發(fā)解算激波管的入射波系及其流場(chǎng)參數(shù)：（1） 入射激波關(guān)系由跨間斷的質(zhì)量、動(dòng)量、能量守恒方程外加PR 方程組成，4 個(gè)方程可建立波后壓強(qiáng)p、比容v、溫度T 和速度u 等4 個(gè)未知數(shù)和入射激波強(qiáng)度的關(guān)系。PR 方程和平衡熱力學(xué)參數(shù)主要通過(guò)能量方程中的焓h 發(fā)生作用。間斷關(guān)系方程全部為代數(shù)方程，因而以牛頓迭代數(shù)值求解（由已知的區(qū)域①參數(shù)出發(fā)，獲得區(qū)域②的參數(shù)）方程組：

式中：下標(biāo)為1 的表示波前氣體參數(shù)，下標(biāo)為2 的表示波后氣體參數(shù)。

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

式中：s、cV 分別為氣體的比熵、比定容熱容，cV，IG 為理想氣體狀態(tài)下氣體的比定壓熱容，c 和u 為當(dāng)?shù)芈曀俸彤?dāng)?shù)厮俣取?/p>

區(qū)域②和區(qū)域③由接觸面隔開，因此區(qū)域②與區(qū)域③的速度（u）和壓強(qiáng)（p）應(yīng)當(dāng)分別相等。可在pu圖上尋找激波管問(wèn)題的解點(diǎn)。

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

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

為呈現(xiàn)廣泛條件下驅(qū)動(dòng)段、被驅(qū)動(dòng)段間壓強(qiáng)配比與入射激波間的關(guān)系，按照以上理論方法對(duì)2 種狀態(tài)方程下不同高、低壓條件驅(qū)動(dòng)產(chǎn)生的入射激波強(qiáng)度進(jìn)行計(jì)算。圖17 所示為以壓強(qiáng)200 和50 MPa 氫氣驅(qū)動(dòng)壓強(qiáng)10～1 000 kPa 的氮?dú)獾娜肷浼げR赫數(shù)，其中空心點(diǎn)為理想氣體狀態(tài)方程的計(jì)算結(jié)果，實(shí)心點(diǎn)為PR 狀態(tài)方程的計(jì)算結(jié)果。由圖17 可以看到：首先，考慮高壓效應(yīng)的入射激波馬赫數(shù)要低于理想氣體情形，這一規(guī)律在全范圍內(nèi)維持不變；其次，高壓真實(shí)氣體效應(yīng)對(duì)入射激波馬赫數(shù)的影響幅度隨驅(qū)動(dòng)段壓強(qiáng)和高、低壓段壓比的增大而增大，但整體變化幅度并不十分顯著，在所測(cè)試的最極端情形下（200 MPa 氫氣驅(qū)動(dòng)10 kPa 氮?dú)猓瑝罕?0 000），理想氣體狀態(tài)方程下入射激波馬赫數(shù)約為11.37，PR 方程下激波馬赫數(shù)約為10.78，兩者差異僅為0.55。

2.3.3 激波管后續(xù)波系傳播

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

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

由于喉道處于壅塞狀態(tài)，反射激波作用后接觸面的運(yùn)動(dòng)速度很大程度上取決于駐室氣流的狀態(tài)。根據(jù)一維等熵流理論，壅塞流速大致正比于駐室區(qū)的聲速（或溫度0.5 次方），由于理想氣體中區(qū)域⑤的氣流溫度略高于真實(shí)氣體的（圖12），因此真實(shí)氣體接觸面的右行速度本來(lái)就略低于理想氣體，但兩者差異并不明顯。實(shí)際流動(dòng)中真正使得真實(shí)氣體中接觸面延遲抵達(dá)的原因是反射稀疏波提前抵達(dá)并作用于駐室區(qū)，使得當(dāng)?shù)芈曀偌眲〗档停瑥亩鴺O大降低了接觸面的后續(xù)運(yùn)動(dòng)速度。

2.4 不同條件下高壓效應(yīng)對(duì)試驗(yàn)時(shí)間的影響

激波風(fēng)洞試驗(yàn)往往需要在同一設(shè)備上獲取不同總溫和不同總壓的試驗(yàn)氣流。通過(guò)改變激波風(fēng)洞運(yùn)行條件，獲得試驗(yàn)氣流總溫2000～6 000 K、總壓20～200 MPa 范圍內(nèi)縫合接觸面運(yùn)行的不同試驗(yàn)工況，以此考察不同條件下高壓真實(shí)氣體效應(yīng)對(duì)試驗(yàn)時(shí)間的影響。試驗(yàn)氣流總溫的調(diào)整和縫合條件的實(shí)現(xiàn)通過(guò)調(diào)整驅(qū)動(dòng)段組分（氫氣與氮?dú)獾呐浔龋┖透摺⒌蛪憾螇罕葘?shí)現(xiàn)，試驗(yàn)氣流總壓的改變通過(guò)在維持高、低壓段壓比情況下調(diào)整壓強(qiáng)實(shí)現(xiàn)。各段初始?xì)怏w溫度恒為300 K。

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

在理想氣體狀態(tài)方程下，由于激波管初始?jí)罕群蜏囟染S持不變，初始入射激波和稀疏波的速度均不變，而在湍流區(qū)激波沿程衰減幅度同樣變化不大，因此，各波系的時(shí)空關(guān)系基本維持不變，反射激波、反射稀疏波和接觸面的抵達(dá)時(shí)刻均無(wú)顯著改變。此時(shí)各工況均為接觸面率先抵達(dá)，有效試驗(yàn)由t2?t1 決定。

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

進(jìn)一步對(duì)試驗(yàn)氣流總溫2 000～6 000 K、總壓20～200 MPa 范圍的其余工況進(jìn)行計(jì)算。為達(dá)到縫合接觸面運(yùn)行條件，不同試驗(yàn)氣流名義總溫對(duì)應(yīng)激波管驅(qū)動(dòng)段摻入氮?dú)饽柗謹(jǐn)?shù)以及高、低壓段的壓比如表4 所示。

將以上測(cè)試范圍內(nèi)高壓真實(shí)氣體效應(yīng)對(duì)有效試驗(yàn)時(shí)間的影響情況總結(jié)為圖20，其中Δ 定義為2 種狀態(tài)方程下有效試驗(yàn)時(shí)間的差值相對(duì)理想氣體方程有效試驗(yàn)時(shí)間的比：

分析固定總壓條件下，有效試驗(yàn)時(shí)間隨總溫的變化情況可知，隨著總溫升高， Δ 逐漸減小。結(jié)合2.2 節(jié)與2.3 節(jié)中的分析，溫度的影響主要體現(xiàn)在試驗(yàn)氣流參數(shù)上，對(duì)于試驗(yàn)時(shí)間影響較小。Δ 隨總溫升高逐漸減小的原因?yàn)楦邷貙?duì)高壓氣體效應(yīng)的影響具有抑制作用。隨著溫度的升高，高壓氣體效應(yīng)引起的稀疏波運(yùn)動(dòng)速度變化幅度減小，在限制有效試驗(yàn)時(shí)間的因素為反射稀疏波到達(dá)后，理想氣體狀態(tài)方程和PR 狀態(tài)方程下有效試驗(yàn)時(shí)間的差異進(jìn)一步減小。

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

綜上所述，對(duì)于高壓氫氣驅(qū)動(dòng)的激波風(fēng)洞內(nèi)的全場(chǎng)流動(dòng)時(shí)空結(jié)構(gòu)和駐室區(qū)氣流參數(shù)的分析獲得如下認(rèn)知：（1） 高壓真實(shí)氣體效應(yīng)對(duì)激波管流動(dòng)的影響主要體現(xiàn)為稀疏波傳播速度加快侵蝕有效試驗(yàn)時(shí)間長(zhǎng)度；（2） 同等壓強(qiáng)下，溫度越高則高壓氣體效應(yīng)越弱。因此，為了減小高壓真實(shí)氣體效應(yīng)的影響，工程中易于實(shí)現(xiàn)的方法有：①延長(zhǎng)驅(qū)動(dòng)段的長(zhǎng)度，使得稀疏波需要經(jīng)過(guò)更長(zhǎng)的距離才能反射到達(dá)駐室，延長(zhǎng)有效試驗(yàn)時(shí)間；②對(duì)激波管驅(qū)動(dòng)段進(jìn)行加熱處理（同氣體介質(zhì)、同等壓強(qiáng)下，高溫氣體驅(qū)動(dòng)能力更強(qiáng)，氣體密度更小，可削弱高壓氣體效應(yīng)）。

3 結(jié) 論

推導(dǎo)了立方型氣體狀態(tài)方程下氣體的熱力學(xué)參數(shù)（內(nèi)能、比熱容、比熱比等）和聲速，并將上述參數(shù)運(yùn)用于耦合高溫?zé)峄瘜W(xué)非平衡過(guò)程的準(zhǔn)一維數(shù)值模擬，對(duì)考慮高壓真實(shí)氣體效應(yīng)的高焓激波風(fēng)洞流動(dòng)結(jié)構(gòu)以及試驗(yàn)氣流參數(shù)進(jìn)行了準(zhǔn)一維數(shù)值計(jì)算和分析，獲得了如下主要結(jié)論。

（1） 采用PR 氣體狀態(tài)方程可以很好地刻畫氣體的高壓真實(shí)氣體效應(yīng)，耦合該方程的準(zhǔn)一維數(shù)值模擬可以準(zhǔn)確復(fù)現(xiàn)高焓激波風(fēng)洞的整體流動(dòng)。

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

（3） 高壓真實(shí)氣體效應(yīng)對(duì)激波風(fēng)洞駐室試驗(yàn)氣流的影響主要呈現(xiàn)為它對(duì)有效試驗(yàn)時(shí)間長(zhǎng)度的侵蝕。這種作用主要通過(guò)高壓效應(yīng)加快稀疏波的傳播、從而致使稀疏波提前抵達(dá)駐室區(qū)并顯著降低當(dāng)?shù)貧饬骺倝哼_(dá)成。如高壓效應(yīng)不能使反射稀疏波先于接觸面抵達(dá)喉道，則高壓效應(yīng)不對(duì)有效試驗(yàn)時(shí)間形成實(shí)質(zhì)侵蝕。

（4） 激波風(fēng)洞在試驗(yàn)氣流總溫較低而總壓較高的工況下，有效試驗(yàn)時(shí)間一般以反射稀疏波頭抵達(dá)時(shí)刻為終點(diǎn)，高壓真實(shí)氣體效應(yīng)更容易在此類工況下對(duì)有效試驗(yàn)時(shí)間產(chǎn)生大幅影響。對(duì)于采用冷高壓氣體驅(qū)動(dòng)的高焓激波風(fēng)洞，可在理想設(shè)計(jì)的基礎(chǔ)上適當(dāng)延長(zhǎng)驅(qū)動(dòng)段長(zhǎng)度以抵消反射稀疏波提前抵達(dá)對(duì)試驗(yàn)時(shí)間的影響。

參考文獻(xiàn)：

[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 風(fēng)洞Ma=10 高超聲速推進(jìn)試驗(yàn)技術(shù)探索 [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]. 宇航學(xué)報(bào)， 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]侯自豪， 朱雨建， 薄靖龍， 等. 真空管道列車準(zhǔn)一維氣動(dòng)特性 [J]. 機(jī)械工程學(xué)報(bào)， 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]中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所， 中國(guó)科學(xué)技術(shù)大學(xué). 超高壓激波類設(shè)備熱化學(xué)非平衡準(zhǔn)一維計(jì)算軟件. 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.

（責(zé)任編輯 張凌云）

基金項(xiàng)目： 國(guó)家自然科學(xué)基金（12172354，U21B6003）