動(dòng)力學(xué)格式在超聲速非平衡流數(shù)值模擬中的應(yīng)用

卓長飛，武曉松，封 鋒

（南京理工大學(xué) 機(jī)械工程學(xué)院，南京210094）

高超聲速飛行器在再入大氣層的過程中與大氣層相互作用，在其頭部周圍形成一道強(qiáng)弓形激波，激波后高溫氣體發(fā)生化學(xué)反應(yīng)，影響到飛行器周圍流場的氣動(dòng)特性和物面熱交換。再入飛行器軌道上的大部分區(qū)域，流場處于化學(xué)非平衡狀態(tài)。因此，對于高超鈍頭體飛行器模擬主要采用有限速率化學(xué)非平衡流模擬。國內(nèi)外在超聲速化學(xué)非平衡流方面的研究已經(jīng)取得一定進(jìn)展，發(fā)展了一些高效、高精度的數(shù)值方法。

在計(jì)算流體力學(xué)中常見的求解方法主要分為兩大類［1］：Godunov方法和Boltzmann方法。Godunov方法是基于黎曼解構(gòu)造通量；Boltzmann方法是基于微觀粒子分布函數(shù)構(gòu)造通量，這種基于氣體分子動(dòng)力學(xué)理論求解可壓縮流動(dòng)的數(shù)值格式稱為GKS（Gas Kinetic Scheme）。由于氣體動(dòng)力學(xué)格式具有健壯性和高效性，在近幾十年已經(jīng)得到廣泛的應(yīng)用，并且在不斷發(fā)展之中。目前主要有兩大類氣體動(dòng)力學(xué)格式：KFVS格式（Kinetic Flux－vector Splitting Scheme）和BGK格式（Bhatnagar Gross Krook Scheme）。KFVS格式是基于求解最簡單的無碰撞的Boltzmann方程，在采用KFVS格式構(gòu)造數(shù)值通量時(shí)需要計(jì)算exp指數(shù)函數(shù)和erf誤差函數(shù)，因此其計(jì)算效率顯得略低。Agarwal和Acheson［2］提出了一種由現(xiàn)有的KFVS格式和WPS格式結(jié)合而成的高效率高精度的格式Kinetic Wave／Particle Splits Scheme（KWPS格式）。從本質(zhì)上來說KWPS格式也由氣體動(dòng)力學(xué)理論推導(dǎo)而來，因此它保留了KFVS格式的精確性和健壯性，同時(shí)在構(gòu)造數(shù)值通量過程中不需要計(jì)算exp指數(shù)函數(shù)和erf誤差函數(shù)，減小了計(jì)算量，保持了WPS格式的高效性。

國內(nèi)外文獻(xiàn)表明氣體動(dòng)力學(xué)格式得到廣泛的發(fā)展和應(yīng)用［1－8］，但用于化學(xué)非平衡流數(shù)值模擬的文獻(xiàn)較少［4］。為了驗(yàn)證2種動(dòng)力學(xué)格式（KFVS格式和KWPS格式）能否在超聲速化學(xué)非平衡流場數(shù)值模擬中仍然表現(xiàn)出高精度、高分辨率特性，同時(shí)驗(yàn)證這2種氣體動(dòng)力學(xué)格式能否應(yīng)用于非結(jié)構(gòu)網(wǎng)格中，本文擬將在結(jié)構(gòu)網(wǎng)格提出的2種氣體動(dòng)力學(xué)格式（KFVS格式和KWPS格式），結(jié)合物理量線性重構(gòu)技術(shù)和限制器，得到基于非結(jié)構(gòu)網(wǎng)格高精度、高分辨率的動(dòng)力學(xué)通量分裂格式有限體積方法，并推廣到超聲速化學(xué)非平衡流場數(shù)值模擬中，為計(jì)算超聲速化學(xué)非平衡流提供一種新方法。

1 數(shù)值方法

1.1 控制方程

在笛卡爾坐標(biāo)系下微分守恒形式的二維軸對稱化學(xué)非平衡流Euler方程為

式中：Q為守恒變量；F，G分別為軸向和徑向無粘通量；S1為軸對稱源項(xiàng)；S0為化學(xué)反應(yīng)源項(xiàng)。這里僅列出主要的變量：

式中：ρ為密度；u，v分別為軸向、徑向速度；p為壓力；E為單位體積的總能；wl（l＝1，…，N）為l組分質(zhì)量分?jǐn)?shù)，ωl為l組分生成速率，N為總組分?jǐn)?shù)。

1.2 動(dòng)力學(xué)格式

本文采用格心格式的有限體積法，因此首先對微分方程進(jìn)行體積分，再應(yīng)用散度定理得到積分形式的控制方程：

式中：Wk，Ak分別為控制體單元第k個(gè)邊界的通量、面積；Nf為控制體單元的邊界個(gè)數(shù)。

計(jì)算控制體單元第k個(gè)邊界對流通量Wk的KFVS格式表達(dá)式為

式中：

計(jì)算控制體單元第k個(gè)邊界對流通量Wk的KWPS格式表達(dá)式如下：

式中：

顯然，這2種動(dòng)力學(xué)矢通量分裂格式不包含自由參數(shù)，而且比常見的FVS類、FDS類和AUSM類格式形式簡單，易于編程。

1.3 物理量重構(gòu)

在計(jì)算單元邊界面通量時(shí)，直接取邊界兩邊單元的物理量，在空間上只有一階精度，數(shù)值耗散較大。為了得到高階精度，需要進(jìn)行空間高階插值。本文利用變量線性重構(gòu)技術(shù)［9］獲得空間二階精度。物理量重構(gòu)表達(dá)式為

式中：U表示守恒變量；i，j分別為左、右單元中心點(diǎn)編號；r為兩點(diǎn)位置矢量；守恒變量的梯度通過格林－高斯公式計(jì)算：

式中：Vi為i單元的體積；nij為i，j單元交界面的單位法矢量；Aij為i，j單元交界面的面積。

1.4 限制器

在跨聲速和超聲速計(jì)算中，為了抑制激波附近可能出現(xiàn)的非物理振蕩，保持?jǐn)?shù)值格式的單調(diào)性，需要使用限制器。限制器通過對物理量的梯度進(jìn)行限制，從而避免在重構(gòu)過程中出現(xiàn)新的極值。本文采用收斂性較好的Venkatarishnan限制器，具體表達(dá)式參見文獻(xiàn)［9］，這里不再列出。使用限制器后，物理量重構(gòu)公式可以寫作：

式中：φi，φj分別表示i，j單元內(nèi)的限制器。

1.5 化學(xué)反應(yīng)模型

化學(xué)反應(yīng)模型是決定成功模擬化學(xué)非平衡流的關(guān)鍵之一。在本文的驗(yàn)證算例中，氫氣和氧氣的化學(xué)反應(yīng)機(jī)理采用7組分8步基元反應(yīng)模型［10］，氮?dú)夂脱鯕獾幕瘜W(xué)反應(yīng)機(jī)理采用高溫空氣5組分12反應(yīng)模型［11］。詳細(xì)反應(yīng)表達(dá)式與參數(shù)如表1和表2所示。表中：阿烏累里斯公式表達(dá)式為Kf＝ATbexp（－E0R0－1T－1）；指 前 因 子A的 單 位為（cm3／mole）m－1·s－1，其中m為基元反應(yīng)級數(shù)，b為溫度指數(shù)，E0為反應(yīng)活化能，R0為氣體常數(shù)，M表示第三碰撞體。

表1 H2＋O2基元反應(yīng)模型

表2 N2＋O2基元反應(yīng)模型

1.6 時(shí)間算子分裂算法

化學(xué)非平衡流控制方程組分為流動(dòng)部分和化學(xué)反應(yīng)部分，兩者相互耦合。本文采用時(shí)間算子分裂算法處理這種耦合過程［12］，即把反應(yīng)流方程組分為流動(dòng)和化學(xué)反應(yīng)2部分，如式（8）和式（9）所示，并把求解流動(dòng)偏微分方程時(shí)采用的時(shí)間步長進(jìn)一步細(xì)分，作為求解化學(xué)反應(yīng)剛性常微分方程的步長，計(jì)算化學(xué)反應(yīng)對流場的貢獻(xiàn)。具體做法是先凍結(jié)化學(xué)反應(yīng)求解得到流場參數(shù)，然后將化學(xué)反應(yīng)看做等容放熱或吸熱過程，保持內(nèi)能、速度參數(shù)不變，計(jì)算各組分的質(zhì)量變化率，最后迭代求解溫度得到流場其他參數(shù)。

2 計(jì)算結(jié)果與分析

2.1 H2＋O2激波誘導(dǎo)燃燒

當(dāng)按理想化學(xué)當(dāng)量比的混合氫氣／氧氣或者氫氣／空氣以超聲速或高超聲速流過鈍頭體時(shí)，會產(chǎn)生較強(qiáng)的弓形激波，波后溫度升高會引起混合氣體燃燒，這就對所采用的數(shù)值計(jì)算方法精度提出了更高要求。Lehr對此進(jìn)行了大量實(shí)驗(yàn)研究，獲得激波形狀、激波位置的數(shù)據(jù)。研究人員又進(jìn)行數(shù)值模擬，得到與實(shí)驗(yàn)比較吻合的結(jié)果。本文選用其中2個(gè)工況驗(yàn)證本文數(shù)值方法的可靠性。

自由來流條件：來流馬赫數(shù)5.08，靜溫291.5K，來流速度2 705.0m／s，靜壓24 797Pa，氫氣和氧氣的體積混合比為2∶1。球頭半徑R＝0.007 5m，參考的實(shí)驗(yàn)數(shù)據(jù)取自文獻(xiàn)［13］，參考的數(shù)值模擬數(shù)據(jù)取自文獻(xiàn)［14－15］。非結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)總數(shù)為12 409，單元邊總數(shù)為36 578，單元總數(shù)為24 170。

圖1給出了2種動(dòng)力學(xué)格式計(jì)算得到的密度等值線。激波與燃燒波耦合形成爆震波，激波層厚度略有增加，2種格式得到的激波面位置均與實(shí)驗(yàn)數(shù)據(jù)［13］基本吻合。圖2～圖4分別給出了2種格式計(jì)算得到的駐點(diǎn)線（無量綱距離X＊：駐點(diǎn)線上距駐點(diǎn)距離與球頭半徑之比）上無量綱密度ρ＊（駐點(diǎn)線上密度與來流密度之比）、無量綱壓力p＊（駐點(diǎn)線上靜壓與來流靜壓之比）與無量綱溫度T＊（駐點(diǎn)線上靜溫與來流靜溫之比）和主要組分質(zhì)量分?jǐn)?shù)w分布，所有計(jì)算結(jié)果與參考文獻(xiàn)［14－15］相比基本吻合。其中，本文計(jì)算得到的密度和壓力與參考文獻(xiàn)計(jì)算的結(jié)果相比，受激波與燃燒波耦合影響而產(chǎn)生的Von Neumann尖峰更為陡峭，說明了本文的數(shù)值模擬更加精確，同時(shí)說明了動(dòng)力學(xué)格式在超聲速化學(xué)非平衡流中高精度、高分辨率特性，精確地分辨了流場中的細(xì)致結(jié)構(gòu)。

圖1 H2＋O2鈍體繞流場密度等值線圖

圖2 H2＋O2鈍體繞流場駐點(diǎn)線上無量綱密度分布

圖3 H2＋O2鈍體繞流場駐點(diǎn)線上無量綱壓力和無量綱溫度分布

圖4 H2＋O2鈍體繞流場駐點(diǎn)線上主要組分質(zhì)量分?jǐn)?shù)分布

2.2 H2＋Air激波誘導(dǎo)燃燒

自由來流條件：來流馬赫數(shù)6.46，靜溫286.6K，來流速度2 605.0m／s，靜壓42 662Pa，氫氣、氧氣、氮?dú)獾捏w積混合比為2∶1∶3.76，在本算例的波后溫度條件下可不計(jì)N2與O2之間的反應(yīng)。球頭半徑R＝0.007 5m，參考的實(shí)驗(yàn)數(shù)據(jù)取自文獻(xiàn)［13］，參考的數(shù)值模擬數(shù)據(jù)取自文獻(xiàn)［14－15］。非結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)總數(shù)為9 663，單元邊總數(shù)為28 440，單元總數(shù)為18 778。

圖5給出了2種動(dòng)力學(xué)格式計(jì)算得到的密度等值線。由于N2的稀釋作用，燃燒釋放出的熱量小于自由來流動(dòng)能，此時(shí)雖然激波與燃燒波重合，但整個(gè)流場中不存在較強(qiáng)的爆震波，2種格式得到的激波面位置均與實(shí)驗(yàn)數(shù)據(jù)［13］基本吻合。圖6、圖7和圖8分別給出了2種格式計(jì)算得到的駐點(diǎn)線上無量綱密度ρ＊、無量綱壓力p＊與無量綱溫度T＊和主要組分質(zhì)量分?jǐn)?shù)w分布。由于N2的稀釋作用，壓力不再出現(xiàn)Von Neumann尖峰，與參考文獻(xiàn)［14－15］的計(jì)算結(jié)果相比，本文計(jì)算得到的Von Neumann尖峰仍然更為陡峭，說明了動(dòng)力學(xué)格式在該算例中精確描述了激波與燃燒波耦合產(chǎn)生的陡峭的Von Neumann尖峰。需要說明的是，國內(nèi)多篇文獻(xiàn)采用不同數(shù)值方法模擬球頭激波燃燒算例均采用文獻(xiàn)［14－15］作為參考，計(jì)算得到的Von Neumann尖峰均比參考文獻(xiàn)陡峭。這可能是參考文獻(xiàn)年代比較久遠(yuǎn)，限于當(dāng)時(shí)數(shù)值模擬的水平，數(shù)值模擬結(jié)果沒能較好地描述Von Neumann尖峰。

圖5 H2＋Air鈍體繞流場密度等值線圖

圖6 H2＋Air鈍體繞流場駐點(diǎn)線上無量綱密度分布

圖7 H2＋Air鈍體繞流場駐點(diǎn)線上無量綱壓力和無量綱溫度分布

圖8 H2＋Air鈍體繞流場駐點(diǎn)線上主要組分質(zhì)量分?jǐn)?shù)分布

2.3 N2＋O2球錐高超聲速非平衡流場

KIM［16］推導(dǎo)了AUSMPW＋格式并用于數(shù)值模擬高超聲速空氣通過球錐的非平衡流場，與實(shí)驗(yàn)結(jié)果對比表明AUSMPW＋格式在高超非平衡流數(shù)值模擬中表現(xiàn)效果較好。本文采用該算例作為驗(yàn)證數(shù)值方法的可靠性算例之一，化學(xué)動(dòng)力學(xué)模型采用高溫空氣5組分12反應(yīng)模型。

自由來流條件：來流速度3 630m／s，靜溫293K，靜壓2 400Pa。球錐半徑R＝0.007m，參考的實(shí)驗(yàn)數(shù)據(jù)取自文獻(xiàn)［16］。非結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)總數(shù)為16 145，單元邊總數(shù)為48 001，單元總數(shù)為31 856。文獻(xiàn)［16］的數(shù)值模擬結(jié)果中沒有給出駐點(diǎn)線上參數(shù)分布，因此本文給出駐點(diǎn)線上主要參數(shù)分布。

圖9中給出了計(jì)算得到的密度、壓力等值線，其中每部分的上半部分為KFVS格式計(jì)算結(jié)果，下半部分為KWPS格式計(jì)算結(jié)果。2種格式計(jì)算得到的激波面位置均與文獻(xiàn)［16］給出的實(shí)驗(yàn)數(shù)據(jù)吻合。此外，本算例的來流馬赫數(shù)約為11.0，與文獻(xiàn)［16］采用AUSMPW＋格式計(jì)算結(jié)果相同，并且不會出現(xiàn)文獻(xiàn)中提及的在高馬赫數(shù)下計(jì)算產(chǎn)生的carbuncle現(xiàn)象。

圖10和圖11分別給出了駐點(diǎn)線上2種格式計(jì)算得到的駐點(diǎn)線上的無量綱溫度和主要組分分布。高速空氣通過球錐產(chǎn)生激波，波后溫度迅速升高，N2和O2開始吸收熱量并發(fā)生離解和置換反應(yīng)，由于N2與O2反應(yīng)過程吸收熱量，因此沿駐點(diǎn)線上溫度逐漸降低，N2和O2組分質(zhì)量分?jǐn)?shù)也均降低，NO組分質(zhì)量分?jǐn)?shù)升高。在該算例溫度下，N2離解較少，N2質(zhì)量分?jǐn)?shù)降低的主要原因是發(fā)生置換反應(yīng)生成NO。而O2質(zhì)量分?jǐn)?shù)降低的主要原因一方面是O2離解，另一方面是生成NO。結(jié)合圖10和圖11可以看出，當(dāng)駐點(diǎn)線軸線溫度低于N2離解溫度時(shí)（本算例條件下離解起始溫度為4 000K左右），N2不再離解，僅發(fā)生置換反應(yīng)，但O2的離解和置換反應(yīng)仍會進(jìn)行。因此，在趨近于壁面處溫度和組分質(zhì)量分?jǐn)?shù)變化均減緩。

圖9 空氣高超聲速鈍體繞流場密度和壓力的等值線圖

圖10 空氣高超聲速鈍體繞流場駐點(diǎn)線上無量綱溫度分布

圖11 空氣高超聲速鈍體繞流場駐點(diǎn)線上主要組分質(zhì)量分?jǐn)?shù)分布

從以上3個(gè)算例看出，KFVS和KWPS格式計(jì)算的結(jié)果基本重合，并且均與參考文獻(xiàn)數(shù)據(jù)基本吻合，說明這2種動(dòng)力學(xué)格式在超聲速化學(xué)非平衡流場數(shù)值模擬中均能較好地分辨出激波和燃燒波等物理現(xiàn)象，流場結(jié)構(gòu)清晰，在強(qiáng)間斷處基本無振蕩，沒有出現(xiàn)非物理解，具有較高精度和分辨率。

3 結(jié)論

本文將2種基于氣體分子動(dòng)力學(xué)理論的動(dòng)力學(xué)格式（KFVS格式和KWPS格式）用于基于非結(jié)構(gòu)網(wǎng)格超聲速化學(xué)非平衡流場數(shù)值模擬中，并且成功模擬了幾個(gè)經(jīng)典的超聲速化學(xué)非平衡流算例。針對本文數(shù)值模擬的算例，可以得到以下結(jié)論：

①將原來在結(jié)構(gòu)網(wǎng)格中提出的2種動(dòng)力學(xué)通量格式（KFVS格式和KWPS格式）推廣應(yīng)用于非結(jié)構(gòu)網(wǎng)格中，數(shù)值結(jié)果表明該方法是可行的。

②將2種動(dòng)力學(xué)格式（KFVS格式和KWPS格式）應(yīng)用到超聲速化學(xué)非平衡流場的數(shù)值模擬中，數(shù)值模擬結(jié)果與前人實(shí)驗(yàn)數(shù)據(jù)和數(shù)值模擬結(jié)果進(jìn)行了比較，說明了本文將2種動(dòng)力學(xué)格式推廣到超聲速化學(xué)非平衡流場數(shù)值模擬中是可行的，并且2種動(dòng)力學(xué)格式在超聲速化學(xué)非平衡流場數(shù)值模擬中均保持高精度、高分辨率特性，為計(jì)算超聲速化學(xué)非平衡流提供一種新方法。

③2種動(dòng)力學(xué)格式不含自由參數(shù)，形式簡單，同時(shí)能夠正確地分辨流場中的激波、燃燒波和其他物理現(xiàn)象，流場結(jié)構(gòu)清晰，具有較高精度和高分辨率；2種動(dòng)力學(xué)格式在超聲速化學(xué)非平衡流場中計(jì)算精度和分辨率基本一致。

［1］LI S H，XU K.Entropy analysis of kinetic flux vector splitting schemes for the compressible Euler equations［J］.Zeitschrift f¨ur Angewandte Mathematik und Physik，2001，52：62－78.

［2］ACHESON K E，AGARWAL R K.A kinetic theory－based wave／particle flux－splitting scheme for the Euler equations.AIAA95－2178［R］.1995.

［3］REKSOPRODJO H S R，AGARWAL R.Implicit kinetic schemes for Euler equations，AIAA2001－2629［R］.2001.

［4］EPPARD W M，GROSSMAN B.An upwind kinetic flux－vector splitting method for flows chemical and thermal non－equilibrium，AIAA93－0894［R］.1993.

［5］RAVICHANDRAN K S.Higher order KFVS algorithms using compact upwind difference operators［J］.Journal of Computational Physics，1997，130：161－173.

［6］MANDAL J C，JAIN K.A new implicit formulation of KFVSscheme for Euler equations，AIAA2006－3 709［R］.2006.

［7］LEI T.Progress in gas－kinetic upwind schemes for the solution of Euler／Navier－Stokes equations－I：overview［J］.Computer ＆Fluids，2012，56：39－48.

［8］湯華中，鄔華謨.高分辨率KFVS有限體積方法及其CFD應(yīng)用［J］.計(jì)算數(shù)學(xué)，1999，21（3）：375－384.TANG Hua－zhong，WU Hua－mo.High resolution KFVS finite volume methods and its appliciton in CFD［J］.Mathematica Numerica Sinica，1999，21（3）：375－384.（in Chinese）

［9］VENKATARISHNAN V.On the accuracy of limiters and convergence to steady state solutions，AIAA93－0880［R］.1993.

［10］劉君，周松柏，徐春光.超聲速流動(dòng)中燃燒現(xiàn)象的數(shù)值模擬方法及應(yīng)用［M］.長沙：國防科技大學(xué)出版社，2008：288－289.LIU Jun，ZHOU Song－bai，XU Chun－guang.The application and method of combustion phenomena in supersonic flow［M］.Changsha：National University of Defense Technology Press，2008：288－289.（in Chinese）

［11］PARK C，YOON S.Caluclation of real－gas effects on blunt－body trim angles［J］.AIAA Joural，1992，30（4）：999－1 007.

［12］胡湘渝，張德良，姜宗林.氣相爆轟基元反應(yīng)模型數(shù)值模擬［J］.空氣動(dòng)力學(xué)報(bào)，2003，21（1）：59－66.HU Xiang－yu，ZHANG De－liang，JIANG Zong－lin.Numerical simulation of gaseous detonation with detailed chemical reaction model［J］.Acta Aerodynamica Sinica，2003，21（1）：59－66.（in Chinese）

［13］LEHR H F.Experiments on shock－induced combustion［J］.Aeronautica Acta，1972，17：589－597.

［14］SOETRISNO M，IMLAY S T.Simulation of the flow field of a ram aceelerator，AIAA91－1915［R］.1991.

［15］YUNGSTER S，EBERHARDT S，BRUCKNER A P.Numerical simulation of hypervelocity projectiles in detonable gases［J］.AIAA Journal，1991，29（2）：187－192.

［16］KIM K H.Accurate computation of hypersonic flows using AUS－MPW ＋ scheme and shock－induced grid technique，AIAA98－2442［R］.1998.