楊天鵬，嚴(yán)東升，于江鵬，郭 陽(yáng)，劉文伶

沿駐點(diǎn)線化學(xué)非平衡程度評(píng)估方法

楊天鵬，嚴(yán)東升，于江鵬，郭 陽(yáng)，劉文伶

（北京航天長(zhǎng)征飛行器研究所，北京，100076）

為了研究高溫真實(shí)氣體效應(yīng)中的化學(xué)非平衡的程度，針對(duì)球頭的沿駐點(diǎn)線流動(dòng)，通過(guò)高溫高壓氣體的化學(xué)反應(yīng)歷程估算化學(xué)反應(yīng)特征時(shí)間，使用正激波流動(dòng)關(guān)系來(lái)估算流動(dòng)特征時(shí)間，建立了化學(xué)非平衡評(píng)估模型，給出了化學(xué)非平衡程度隨馬赫數(shù)和高度的變化關(guān)系。計(jì)算結(jié)果表明，馬赫數(shù)越大、高度越低，化學(xué)反應(yīng)特征時(shí)間比流動(dòng)特征時(shí)間越小，定義的化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間量級(jí)相當(dāng)時(shí)為化學(xué)非平衡流動(dòng)，時(shí)間比大于10則接近化學(xué)凍結(jié)流，時(shí)間比小于0.1則接近化學(xué)平衡流。在設(shè)定的球頭模型下，化學(xué)反應(yīng)平衡流在高度0km、10km、20km、30km、40km分別需要馬赫數(shù)大于8.0、9.5、10.8、12.7、20.0，化學(xué)反應(yīng)凍結(jié)流在高度30km、40km、50km、60km分別需要馬赫數(shù)小于8.4、10.2、12.9、17.5，中間來(lái)流狀態(tài)對(duì)應(yīng)化學(xué)反應(yīng)非平衡狀態(tài)。

數(shù)值模擬；高溫真實(shí)氣體效應(yīng)；化學(xué)非平衡

0 引 言

飛行器在大氣中高速運(yùn)動(dòng)時(shí)，其前緣附近產(chǎn)生的強(qiáng)激波會(huì)使得激波層內(nèi)的空氣處于高溫高壓狀態(tài)，此時(shí)空氣會(huì)發(fā)生一系列復(fù)雜的化學(xué)反應(yīng)和能量變化，這種現(xiàn)象被稱為高溫真實(shí)氣體效應(yīng)。對(duì)于高溫真實(shí)氣體繞流，一般可將流動(dòng)分為3種情況[1]：第1種是凍結(jié)流，此時(shí)流場(chǎng)中的化學(xué)反應(yīng)速率為零，氣體組元的組成不隨空間和時(shí)間變化，但粘性流動(dòng)會(huì)使得氣體組元由于擴(kuò)散而改變；第2種是平衡流，此時(shí)流場(chǎng)中的氣體處于熱力學(xué)平衡狀態(tài)，且流體處于化學(xué)平衡狀態(tài)，此時(shí)流場(chǎng)中的化學(xué)反應(yīng)速率為無(wú)窮大；第3種是非平衡流，流場(chǎng)中氣體的熱力學(xué)性質(zhì)和化學(xué)反應(yīng)狀態(tài)隨時(shí)間空間發(fā)生變化，此時(shí)化學(xué)反應(yīng)或能量變化的特征時(shí)間與流動(dòng)特征時(shí)間處于相近量級(jí)。

對(duì)于高溫真實(shí)氣體效應(yīng)中涉及到的化學(xué)平衡態(tài)問(wèn)題，在多個(gè)文獻(xiàn)中都給出了大致的參考[2,3]，但是這些文獻(xiàn)并沒(méi)有對(duì)平衡態(tài)問(wèn)題進(jìn)行定量分析。文獻(xiàn)[4]認(rèn)為可以通過(guò)達(dá)姆科勒數(shù)am，用于描述同一系統(tǒng)中化學(xué)反應(yīng)相比其他現(xiàn)象的相對(duì)時(shí)間尺度，使用粒子碰撞特征時(shí)間與流場(chǎng)特征時(shí)間尺度之比來(lái)定義流動(dòng)的化學(xué)狀態(tài)。但是該定義是一種寬泛的定性認(rèn)識(shí)，一般認(rèn)為特征時(shí)間之比高于4個(gè)量級(jí)后才能被劃分到另一流態(tài)，不便于工程應(yīng)用。由于化學(xué)反應(yīng)與組元構(gòu)成、能量輸運(yùn)、化學(xué)反應(yīng)機(jī)理等等密切相關(guān)，并且流動(dòng)中還需要關(guān)注湍流、邊界層、壁面催化、燒蝕等一系列與化學(xué)反應(yīng)相耦合的問(wèn)題，其表現(xiàn)出強(qiáng)烈的非線性，建模難度非常高，難以給出明確的公式對(duì)其進(jìn)行計(jì)算。

為了給出高溫真實(shí)氣體效應(yīng)中對(duì)化學(xué)非平衡程度的評(píng)估，本文從化學(xué)反應(yīng)的宏觀表現(xiàn)出發(fā)，使用高溫高壓空氣的化學(xué)反應(yīng)歷程時(shí)間與流動(dòng)特征時(shí)間作為比較，建立化學(xué)非平衡評(píng)估模型，給出化學(xué)非平衡程度隨高度和馬赫數(shù)的變化規(guī)律，為流場(chǎng)化學(xué)反應(yīng)模型選取、非平衡流場(chǎng)計(jì)算結(jié)果不確定度等研究提供參考，提升高速飛行器氣動(dòng)特性預(yù)示能力。

1 沿駐點(diǎn)線化學(xué)反應(yīng)流動(dòng)特征分析

為了定量地建立空氣化學(xué)反應(yīng)非平衡程度評(píng)估模型，需要首先對(duì)相應(yīng)流場(chǎng)的物理特征進(jìn)行分析。本文選取化學(xué)非平衡程度最為嚴(yán)重的駐點(diǎn)線區(qū)域作為評(píng)估對(duì)象，使用數(shù)值模擬方法對(duì)化學(xué)反應(yīng)效應(yīng)和流動(dòng)特征進(jìn)行分析。

本文使用多組元可壓縮雷諾平均Navier-Stokes方程作為控制方程，其積分形式[5]為

計(jì)算模型選用半徑50 mm的球頭，網(wǎng)格量為100×100，計(jì)算網(wǎng)格如圖1所示。邊界設(shè)置上來(lái)流為遠(yuǎn)場(chǎng)邊界，物面使用無(wú)滑移、絕熱壁和非催化壁邊界條件，排除物面模型對(duì)流動(dòng)非平衡程度的影響。來(lái)流組元為空氣，使用Gupta給出的5組元六反應(yīng)模型及其化學(xué)反應(yīng)速率系數(shù)[6]，化學(xué)反應(yīng)模型如表1所示。

圖1 計(jì)算網(wǎng)格

表1 化學(xué)反應(yīng)方程

Tab.1 Chemical Reaction Equation

序號(hào)化學(xué)反應(yīng)方程式 1O2+M=2O+M 2N2+M=2N+M 3NO+M=N+O+M 4N2+O=NO+N 5NO+O=O2+N 6N2+N=2N+N

下面給出高度30 km條件下，分別為10、15、20的3種情況下沿駐點(diǎn)線化學(xué)反應(yīng)能量釋放情況。這里定義反應(yīng)的化學(xué)反應(yīng)功率密度E為

式中 fi為化學(xué)反應(yīng)i的凈反應(yīng)速率；ΔHi為化學(xué)反應(yīng)i的焓變。通過(guò)該表達(dá)式可以得知在空間中某化學(xué)反應(yīng)所造成的產(chǎn)熱功率，其反映了某化學(xué)反應(yīng)對(duì)流場(chǎng)內(nèi)能的影響程度，正值為吸熱，負(fù)值為放熱。從圖2中可以看出，不同來(lái)流馬赫數(shù)下，化學(xué)反應(yīng)流場(chǎng)表現(xiàn)出不同的能量釋放模式，在Ma=10的情況下，化學(xué)反應(yīng)吸熱主要以O(shè)2分解為主；在Ma=15的情況下，開(kāi)始出現(xiàn)NO分解；在Ma=20的情況下，NO分解與O2分解并重，且N2分解反應(yīng)也表現(xiàn)出較高的吸熱量。這說(shuō)明在評(píng)估化學(xué)非平衡時(shí)，需要同時(shí)將O2、N2和NO 3個(gè)組元選為評(píng)估指標(biāo)。

續(xù)圖2

圖3給出了=15情況下沿駐點(diǎn)線壓強(qiáng)、溫度和速度分布。可以看到，沿駐點(diǎn)線的壓強(qiáng)在激波后略有升高，基本上保持不變。激波后溫度由于激波層較窄，首先呈現(xiàn)出化學(xué)凍結(jié)的高溫，之后由于空氣分解吸熱導(dǎo)致溫度逐漸降低，最后形成一個(gè)變化較小的溫度平臺(tái)。而激波后的速度則呈現(xiàn)出接近直線的分布。激波后流動(dòng)參數(shù)分布的規(guī)律性使得后文中建立的化學(xué)反應(yīng)非平衡評(píng)估模型更加接近于物理實(shí)際情況。

圖3 Ma=15算例沿駐點(diǎn)線流動(dòng)參數(shù)分布

續(xù)圖3

2 化學(xué)反應(yīng)特征時(shí)間評(píng)估

根據(jù)上文分析得到的沿駐點(diǎn)線流動(dòng)特征，本文使用密閉容器中的化學(xué)反應(yīng)進(jìn)程來(lái)評(píng)估化學(xué)反應(yīng)特征時(shí)間。容器內(nèi)的起始溫度和壓強(qiáng)取來(lái)流總溫和總壓，隨著反應(yīng)進(jìn)程的進(jìn)行，容器內(nèi)的溫度逐漸降低，該模型與圖3給出的沿駐點(diǎn)線上流動(dòng)參數(shù)變化過(guò)程較為接近。

圖4給出了取=15、高度30 km條件下的組元變化歷程，給出了以O(shè)2、N2和NO為參考，組元濃度分別達(dá)到80%、90%和99%平衡狀態(tài)情況下的時(shí)刻位置?？梢钥吹綄⒒瘜W(xué)反應(yīng)特征時(shí)間定義為容器中的氧氣和氮?dú)膺_(dá)到80%平衡狀態(tài)時(shí)的反應(yīng)時(shí)間較為合適，之后化學(xué)反應(yīng)處于反應(yīng)緩慢的平臺(tái)區(qū)，化學(xué)反應(yīng)導(dǎo)致的組元變化和能量變化主要集中于前80%部分，若選擇更高的非平衡標(biāo)準(zhǔn)會(huì)使得化學(xué)反應(yīng)特征時(shí)間成倍增加，化學(xué)非平衡評(píng)估模型的敏感性降低。

圖4 Ma=15算例沿駐點(diǎn)線組元變化歷程

圖5給出了化學(xué)反應(yīng)特征時(shí)間隨馬赫數(shù)和高度的變化關(guān)系，可以看到，馬赫數(shù)越大、高度越低，化學(xué)反應(yīng)特征時(shí)間越小，這主要是由于壓強(qiáng)和溫度增加有助于分子碰撞，從而有利于化學(xué)反應(yīng)進(jìn)行。但是壓強(qiáng)增加并不總有利于空氣化學(xué)反應(yīng)進(jìn)行，可以看到=12附近的曲線會(huì)出現(xiàn)轉(zhuǎn)折，在轉(zhuǎn)折兩側(cè)的曲線斜率發(fā)生了明顯變化。這是由于壓強(qiáng)增加會(huì)增大空氣化學(xué)反應(yīng)中各組元的濃度，使得分解反應(yīng)難以進(jìn)行，出現(xiàn)局部溫度壓強(qiáng)增大并不促進(jìn)化學(xué)反應(yīng)發(fā)生的現(xiàn)象。

圖5 化學(xué)反應(yīng)特征時(shí)間

3 流動(dòng)特征時(shí)間評(píng)估

本文以流體微團(tuán)沿駐點(diǎn)線從正激波后運(yùn)動(dòng)到物面的時(shí)間為流動(dòng)特征時(shí)間。首先對(duì)激波脫體距離進(jìn)行估算，球頭激波脫體距離近似公式[7]如下：

圖6給出了不同高度下流動(dòng)特征時(shí)間隨馬赫數(shù)變化情況，由圖6可知，在本文所使用的模型中，馬赫數(shù)對(duì)流動(dòng)特征時(shí)間起主導(dǎo)作用，不同高度的流動(dòng)特征時(shí)間變化不大，其整體上表現(xiàn)出隨著馬赫數(shù)增大而減小的趨勢(shì)。

圖6 流動(dòng)特征時(shí)間

4 化學(xué)非平衡程度模型

本文將化學(xué)非平衡定義為化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間處于相近量級(jí)，圖7給出了化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間之比隨馬赫數(shù)的變化?？梢钥吹?，由于不同高度下的流動(dòng)特征時(shí)間較為接近，因此化學(xué)非平衡程度整體變化情況與化學(xué)特征時(shí)間模型基本一致。對(duì)于不同高度，高度越高，流動(dòng)越接近于化學(xué)凍結(jié)流，這是由于低壓使得化學(xué)反應(yīng)更難發(fā)生，化學(xué)反應(yīng)時(shí)間更長(zhǎng)；而對(duì)于低高度，流動(dòng)越接近于化學(xué)平衡流，這是由于高壓促進(jìn)分子碰撞，使得化學(xué)反應(yīng)時(shí)間更短。

圖7 化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間之比

下面進(jìn)一步對(duì)比實(shí)際情況與建模之間的不同。為方便起見(jiàn)，選擇高度40 km、馬赫數(shù)為10、13、20的來(lái)流參數(shù)進(jìn)行計(jì)算，其化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間之比分別為12.3、1.13和0.094。圖8給出了沿駐點(diǎn)線不同組元摩爾分?jǐn)?shù)分布的對(duì)比，其中0D代表化學(xué)反應(yīng)特征時(shí)間的評(píng)估結(jié)果，其零時(shí)間對(duì)應(yīng)于激波面，右側(cè)終點(diǎn)為達(dá)到80%平衡的時(shí)間。從圖8中可以看到，對(duì)于=10的情況，流動(dòng)接近于化學(xué)凍結(jié)流，組元整體變化程度較??；對(duì)于=13的情況，流動(dòng)處于明顯的化學(xué)平衡態(tài)，且在接近物面處組元的分布恰好接近于平衡態(tài)，可以認(rèn)為此時(shí)化學(xué)反應(yīng)時(shí)間與流動(dòng)時(shí)間基本相當(dāng)；對(duì)于=20的情況，可以看到沿駐點(diǎn)線組元很快達(dá)到了平衡態(tài)，流動(dòng)接近于化學(xué)平衡流。由于本文采用了正激波模型，因此數(shù)值模擬得到的結(jié)果相比于建模結(jié)果，其激波后溫度較低，化學(xué)反應(yīng)速率較低。

圖8 沿駐點(diǎn)線組元變化歷程與非平衡建模結(jié)果對(duì)比

續(xù)圖8

以上3個(gè)算例對(duì)比說(shuō)明本文所建立的化學(xué)非平衡程度評(píng)估模型是可行的，可以用于指導(dǎo)化學(xué)非平衡數(shù)值模擬。對(duì)于流動(dòng)中化學(xué)非平衡的判據(jù)則需要根據(jù)數(shù)值模擬的需要進(jìn)行，從結(jié)果來(lái)看，化學(xué)反應(yīng)特征時(shí)間比流動(dòng)特征時(shí)間的值取為0.1～10之間可以認(rèn)為是化學(xué)非平衡，大于10則接近化學(xué)平衡流，小于0.1則接近化學(xué)凍結(jié)流。需要注意的是，在低馬赫數(shù)下化學(xué)非平衡程度呈現(xiàn)出向化學(xué)平衡流移動(dòng)的趨勢(shì)，但是此時(shí)流場(chǎng)中的溫度不夠高，高溫真實(shí)氣體效應(yīng)不明顯，即使是化學(xué)非平衡流也可使用化學(xué)凍結(jié)流模型進(jìn)行計(jì)算。

由于前緣半徑與激波脫體距離近似成正比關(guān)系，前緣半徑越大沿駐點(diǎn)線的流動(dòng)特征時(shí)間越長(zhǎng)，但化學(xué)反應(yīng)特征時(shí)間是不具有幾何相似性的，因此前緣半徑越大，化學(xué)反應(yīng)特征時(shí)間比流動(dòng)特征時(shí)間越小。實(shí)際飛行器前緣半徑的取值與本文所選數(shù)值不會(huì)有量級(jí)上的差異，相比于化學(xué)反應(yīng)特征時(shí)間約3～4個(gè)量級(jí)的跨度，前緣半徑改變帶來(lái)的時(shí)間變化相對(duì)較小。

5 結(jié) 論

本文針對(duì)球頭繞流中的駐點(diǎn)線流動(dòng)，通過(guò)高溫高壓氣體的化學(xué)反應(yīng)歷程給出了化學(xué)反應(yīng)特征時(shí)間，使用正激波流動(dòng)關(guān)系來(lái)估算流動(dòng)特征時(shí)間，建立了化學(xué)非平衡評(píng)估模型，通過(guò)化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間之比給出了流場(chǎng)中的化學(xué)非平衡程度隨馬赫數(shù)和高度的變化結(jié)果，并通過(guò)數(shù)值模擬結(jié)果驗(yàn)證了非平衡評(píng)估結(jié)果的正確性?？捎糜谥笇?dǎo)高速流動(dòng)中的化學(xué)反應(yīng)模型選取。

計(jì)算結(jié)果表明，化學(xué)反應(yīng)特征時(shí)間隨馬赫數(shù)增大而減小，隨高度增加而增加；流動(dòng)特征時(shí)間隨高度變化不大，隨馬赫數(shù)增大而減??；化學(xué)反應(yīng)特征時(shí)間與流動(dòng)特征時(shí)間之比與化學(xué)反應(yīng)特征時(shí)間規(guī)律相近，其隨馬赫數(shù)增大而減小，隨高度增加而增加。與數(shù)值模擬結(jié)果對(duì)比表明，時(shí)間比處于0.1～10之間可以認(rèn)為是化學(xué)非平衡，大于10則接近化學(xué)凍結(jié)流，小于0.1則接近化學(xué)平衡流。在本文設(shè)定的球頭模型下，化學(xué)反應(yīng)平衡流在高度0 km、10 km、20 km、30 km、40 km分別需要馬赫數(shù)大于8、9.5、10.8、12.7、20.0，化學(xué)反應(yīng)凍結(jié)流在高度30 km、40 km、50 km、60 km分別需要馬赫數(shù)小于8.4、10.2、12.9、17.5，中間來(lái)流狀態(tài)對(duì)應(yīng)化學(xué)反應(yīng)非平衡狀態(tài)。

[1] 歐陽(yáng)水吾, 謝中強(qiáng). 高溫非平衡空氣繞流[M]. 北京: 國(guó)防工業(yè)出版社, 2001.

Ouyang Shuiwu, Xie Zhongqiang. High temperature non-equilibrium flow[M]. Beijing: National Defence Industry Press, 2001.

[2] Miller J H. Computational aerothermodynamic datasets for hypersonic heat transfer on reentry vehicles[C]. San Francisco: AIAA Atmospheric Flight Mechanics Conference and Exhibit, 2005.

[3] Tauber M E, Meneses G P. Aerothermodynamics of transatmospheric vehicles[C]. Boston: 4th Thermophysics and Heat Transfer Conference, 1986.

[4] 中國(guó)科學(xué)院. 中國(guó)學(xué)科發(fā)展戰(zhàn)略新型飛行器中的關(guān)鍵力學(xué)問(wèn)題[M]. 北京: 科學(xué)出版社, 2018.

Chinese Academy of Sciences. Chinese discipline development strategy key mechanical problems in new aircraft[M]. Beijing: Science Press, 2018.

[5] 王江峰, 等. 高超聲速?gòu)?fù)雜氣動(dòng)問(wèn)題數(shù)值方法研究進(jìn)展[J]. 航空學(xué)報(bào), 2015, 36(1): 159-175.

Wang Jiangfeng, et al. Progress in numerical simulation techniques of hypersonic aerodynamic problems[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 159-175.

[6] Gupta R N, et al. A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000K[R]. NASA-RP-1232, 1989.

[7] Anderson J D. 高超聲速和高溫氣體動(dòng)力學(xué)（第2版）[M]. 北京: 航空工業(yè)出版社, 2013.

Anderson J D. Hypersonic and high-temperature gas dynamics(Second Edition)[M]. Beijing: Aviation Industry Press, 2013.

Evaluation Method of Chemical Non-equilibrium along Stagnation Line

Yang Tian-peng, Yan Dong-sheng, Yu Jiang-peng, Guo Yang, Liu Wen-ling

(Beijng Institute of Space Long March Vehicle, Beijing, 100076)

In order to study the degree of chemical non-equilibrium in high-temperature real gas effect, this paper begins with the flow of the ball head along the stagnation line to establish a chemical non-equilibrium evaluation model. In this model, the chemical reaction characteristic time is calculated by using the chemical reaction process, and the flow characteristic time is evaluated by using the relationship before and after the normal shock wave. The degree of chemical non-equilibrium variation law with Mach number and height is given. Calculation results show that the larger the Mach number and the lower the height makes the ratio of chemical reaction characteristic time and flow characteristic time smaller. The flow is chemical non-equilibrium when the chemical reaction characteristic time and flow characteristic time is in the near order of magnitude. If the time ratio is greater than 10, the flow is close to the chemical equilibrium flow, and if the time ratio is less than 0.1, the flow is close to the chemical frozen flow. Under this ball head model, the chemical equilibrium flow requires Mach numbers greater than 8.0, 9.5, 10.8, 12.7, and 20.0 at heights of 0 km, 10 km, 20 km, 30 km, and 40 km respectively. The chemical frozen flow requires Mach numbers less than 8.4, 10.2, 12.9, and 17.5 at heights of 30 km, 40 km, 50 km, and 60 km respectively. The intermediate inflow state is chemical non-equilibrium flow.

numerical simulation; high-temperature real gas effect; chemical non-equilibrium

2097-1974(2023)01-0122-05

10.7654/j.issn.2097-1974.20230124

V211.3

A

2022-11-09；

2022-12-04

楊天鵬（1993-），男，博士，工程師，主要研究方向?yàn)轱w行器氣動(dòng)設(shè)計(jì)。

嚴(yán)東升（1963-），男，研究員，主要研究方向?yàn)轱w行力學(xué)設(shè)計(jì)。

于江鵬（1995-），男，博士，工程師，主要研究方向?yàn)轱w行器氣動(dòng)設(shè)計(jì)。

郭 陽(yáng)（1991-），男，博士，工程師，主要研究方向?yàn)轱w行器氣動(dòng)設(shè)計(jì)。

劉文伶（1963-），女，研究員，主要研究方向?yàn)轱w行器總體設(shè)計(jì)。