雙錐/雙楔流動(dòng)中的高溫氣體效應(yīng)仿真模擬

袁軍婭，任 翔，蔡國飆，王偉宗，賀碧蛟

(北京航空航天大學(xué)宇航學(xué)院，北京 100191)

引 言

在大氣層內(nèi)或跨大氣層高超聲速飛行的飛行器，與周圍空氣劇烈摩擦并對(duì)前方空氣強(qiáng)烈壓縮，由于黏性耗散效應(yīng)和激波強(qiáng)烈壓縮，氣體損失的巨大動(dòng)能轉(zhuǎn)變?yōu)閮?nèi)能，造成流場(chǎng)溫度急劇升高，形成高溫高壓的氣動(dòng)熱環(huán)境，引起空氣分子一系列復(fù)雜物理化學(xué)反應(yīng)現(xiàn)象如離解、電離和內(nèi)部能態(tài)激發(fā)等，從而造成高溫真實(shí)氣體效應(yīng)[1]。

在高溫真實(shí)氣體描述方面，計(jì)算流體力學(xué)方法和直接模擬Monte Carlo(direct simulation Monte Carlo，DSMC)方法都得到適當(dāng)發(fā)展并在科學(xué)研究和工程實(shí)踐中得到了廣泛應(yīng)用[2-3]。CFD和DSMC方法對(duì)高溫真實(shí)氣體效應(yīng)的模型化包括熱力學(xué)參數(shù)、輸運(yùn)參數(shù)、熱力學(xué)非平衡過程、化學(xué)非平衡過程的描述。目前廣泛使用的各類模型主要是從常溫向高溫?cái)U(kuò)展、平衡態(tài)向非平衡態(tài)擴(kuò)展，在常溫和平衡態(tài)下相互符合的不同模型在高溫非平衡態(tài)下會(huì)產(chǎn)生明顯的差異，這些差異會(huì)影響流場(chǎng)和氣動(dòng)參數(shù)的模擬預(yù)測(cè)。

目前熱化學(xué)非平衡流CFD模擬，自20世紀(jì)60年代提出Lighthill離解模型、Landau-Teller振動(dòng)松弛模型[4]，到80年代的Park雙溫反應(yīng)模型[5]，這些物理模型由于某些局限性或者缺陷而不斷推舊出新。Park雙溫模型主要關(guān)注振動(dòng)態(tài)的松弛過程，并假設(shè)轉(zhuǎn)動(dòng)態(tài)與平動(dòng)態(tài)時(shí)刻達(dá)到平衡，從而使用兩個(gè)溫度來描述流場(chǎng)的熱運(yùn)動(dòng)狀態(tài)。同時(shí)CFD中大氣組分的化學(xué)反應(yīng)速率包含Dunn等[6]，Gupta等[7]和Park[8]機(jī)理數(shù)據(jù)。

DSMC方法主要使用L-B(Larsen-Borgnakke)[9]和量子L-B模型分別模擬分子轉(zhuǎn)動(dòng)和振動(dòng)態(tài)的激發(fā)與松弛，并先后發(fā)展了總碰撞能量(total collision energy，TCE)模型[10]、振動(dòng)控制離解(vibrationally favored dissociation，VFD)模型[11]、廣義碰撞能量(generalized collision energy，GCE)模型[12]和量子動(dòng)力學(xué)(quantum-kinetic，Q-K)反應(yīng)模型[13-15]來描述化學(xué)反應(yīng)過程。其中Q-K模型不同于其他依賴于宏觀物理量的化學(xué)反應(yīng)模型，是僅基于兩個(gè)碰撞粒子的基本屬性的分子級(jí)化學(xué)模型。Q-K模型對(duì)于分子振動(dòng)能采用簡(jiǎn)單諧振子模型，并綜合考慮分子的離解與振動(dòng)能的激發(fā)。Liechty[16]對(duì)Q-K模型進(jìn)一步發(fā)展并囊括了電子能級(jí)和電離反應(yīng)。

雙錐和雙楔模型是驗(yàn)證激波-激波和激波-邊界層干擾、流動(dòng)分離與渦旋等現(xiàn)象的主要模型，目前學(xué)者對(duì)該模型進(jìn)行了廣泛的理論分析、風(fēng)洞測(cè)量和數(shù)值模擬研究[17-25]。1999年，Olejniczak等[18]采用帶化學(xué)反應(yīng)的N2流進(jìn)行雙楔實(shí)驗(yàn)，并采用量熱氣體模型進(jìn)行數(shù)值仿真，模擬的分離區(qū)比實(shí)驗(yàn)值小。2003年，Nompelis等[26]對(duì)高超聲速雙錐實(shí)驗(yàn)的數(shù)值模擬將模型的傳熱速率高估了約20%。2007年，Hash等[27]研究了化學(xué)反應(yīng)速率對(duì)CUBRC LENS設(shè)備中高焓空氣實(shí)驗(yàn)的影響。2012年，Holden等[22]介紹了在LENS超高速地面試驗(yàn)裝置中為獲得基本測(cè)量數(shù)據(jù)而進(jìn)行的實(shí)驗(yàn)方案，并以此來評(píng)價(jià)和發(fā)展用于預(yù)測(cè)非平衡流動(dòng)化學(xué)、邊界層過渡和激波/湍流邊界層相互作用的模型。2017年，Knight等[28]采用Mach數(shù)為7.1的雙楔形模型，在2.1 MJ/kg和8 MJ/kg 的空氣和氮?dú)庵校瑢?duì)高超聲速激波層流邊界層相互作用的CFD預(yù)測(cè)能力進(jìn)行了評(píng)估。

本文通過可描述內(nèi)能松弛和化學(xué)反應(yīng)的CFD和DSMC對(duì)3個(gè)具有試驗(yàn)數(shù)據(jù)的雙錐雙楔流動(dòng)進(jìn)行仿真分析，研究高溫真實(shí)氣體效應(yīng)對(duì)流動(dòng)預(yù)測(cè)的影響。

1 物理模型和數(shù)值方法

1.1 物理模型

使用CFD和DSMC方法分別對(duì)焓值為2.23，8.17，26.1 MJ/kg這3組雙錐或雙楔流動(dòng)進(jìn)行仿真[18，21，28]，來流工質(zhì)均為N2，3組幾何模型見圖1，來流參數(shù)見表1。仿真過程中考慮稀薄效應(yīng)和高溫真實(shí)氣體效應(yīng)的影響。

表1 來流參數(shù)和壁面條件Table 1 Incoming flow parameters and wall conditions

(a) Double cone

1.2 數(shù)值方法

1.2.1 CFD方法

采用Park雙溫模型來模擬高超聲速流動(dòng)中的熱化學(xué)非平衡效應(yīng)，控制方程除包含連續(xù)性方程、動(dòng)量方程、總能量方程外，還添加了振動(dòng)能守恒方程

式中，U為守恒量，通量分為無黏通量Fi，inv和有黏通量Fi，vis,W為源項(xiàng)，下標(biāo)i指代三維空間，對(duì)于Park雙溫模型

U=(ρ，ρs，ρu1，ρu2，ρu3，Eve，E)
W=(0，ωs，0，0，0，ωv，0)

式中，ρ為密度;u為速度;E和Eve分別為總能量和振動(dòng)-電子能;ωs為組分s的質(zhì)量變化，取決于化學(xué)反應(yīng)過程;ωv為振動(dòng)能的變化量，由振動(dòng)-平動(dòng)松弛過程和化學(xué)反應(yīng)共同決定

ωv=QV-T+QC-V

式中，振動(dòng)-平動(dòng)項(xiàng)QV-T采用Landau-Teller公式計(jì)算[4]

振動(dòng)-平動(dòng)松弛時(shí)間τV-T取為Millikan-White擬合與Park修正項(xiàng)的和值[29-30]。

化學(xué)反應(yīng)帶來的振動(dòng)-電子能量變化為

其中，A,β,Ta分別為指前因子、溫度指數(shù)、活化溫度，Tc，f為正向反應(yīng)的控制溫度，由Ttr和Tve共同決定。

這里只考慮N2的離解反應(yīng)(N2+ N2= 2N + N2，N2+ N = 2N + N)且其反應(yīng)速率取自DSMC的Q-K模型[31]，忽略電離反應(yīng)。使用基于第一激發(fā)能級(jí)的簡(jiǎn)單諧振子模型和多電子激發(fā)態(tài)模型共同來確定氣體的熱力學(xué)參數(shù)。

1.2.2 DSMC方法

DSMC方法是由Bird提出的基于概率的直接分子模擬，并將分子運(yùn)動(dòng)與碰撞過程解耦處理[32]。本文使用量子L-B模型模擬分子振動(dòng)態(tài)的激發(fā)與松弛過程。DSMC分別模擬轉(zhuǎn)動(dòng)能和振動(dòng)能的松弛，而CFD的雙溫模型中假設(shè)分子的平動(dòng)和轉(zhuǎn)動(dòng)能間是平衡的。

采用Q-K模型模擬N2的離解反應(yīng)[14]，即對(duì)于一對(duì)碰撞分子，由碰撞能量確定可達(dá)到的最大振動(dòng)能級(jí)

如果此能級(jí)高于分解能級(jí)，則發(fā)生離解。

1.2.3 模型設(shè)置

本文分別在CFD和DSMC方法中考慮高溫條件下振動(dòng)能和化學(xué)反應(yīng)過程的模擬，設(shè)置了如表2中所示的幾種模式，并且為了在結(jié)果中使用方便，約定了簡(jiǎn)寫標(biāo)記。受限于存儲(chǔ)容量和計(jì)算能力，僅對(duì)低焓流動(dòng)使用DSMC方法。

表2 數(shù)值方法設(shè)置Table 2 Numerical method settings

1.3 網(wǎng)格無關(guān)性驗(yàn)證

CFD和DSMC均采用貼體四邊形網(wǎng)格，依據(jù)文獻(xiàn)[33]中的Recell，w=1準(zhǔn)則分別設(shè)定雙錐/雙楔靠近壁面附近的法向網(wǎng)格尺寸。CFD方法具有成熟的離散方法和數(shù)值格式，而DSMC方法對(duì)網(wǎng)格的依賴嚴(yán)格，這里專門對(duì)DSMC方法進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證。

分別采用網(wǎng)格數(shù)量為7.8×104，3.12×105，1.248×106且貼壁網(wǎng)格尺寸為0.5，0.2，0.2 mm的3套網(wǎng)格對(duì)雙錐工況進(jìn)行無關(guān)性驗(yàn)證。圖2展示了獲得的壁面熱流分布。結(jié)果表明，幾套網(wǎng)格主要在分離區(qū)和激波-邊界層干擾位置處存在一定差異，并且后兩套網(wǎng)格的結(jié)果相近，因此在后續(xù)仿真中采用第2套網(wǎng)格。另外兩個(gè)雙楔流動(dòng)也進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證，這里不再贅述。

圖2 網(wǎng)格無關(guān)性Fig.2 Grid independence

2 結(jié)果與討論

2.1 低焓雙錐流動(dòng)

美國 Holden 等的CUBRC第7次試車實(shí)驗(yàn)數(shù)據(jù)作為稀薄流仿真軟件的驗(yàn)證算例。實(shí)驗(yàn)測(cè)量稀薄超聲速來流對(duì)雙圓錐體的氣動(dòng)力和氣動(dòng)熱值，是國際上一個(gè)典型的稀薄氣體動(dòng)力學(xué)實(shí)驗(yàn)。

圖3展示了使用考慮振動(dòng)能松弛及氮?dú)怆x解反應(yīng)的CFD方法對(duì)雙錐流場(chǎng)的模擬結(jié)果。從第一錐的頭尖部產(chǎn)生一道附體激波，因第二錐而形成非附體激波，且形成三波點(diǎn);透射激波入射到第二錐的錐面上，使流動(dòng)分離，且產(chǎn)生一個(gè)很大的當(dāng)?shù)貕荷?，這個(gè)壓升使上游流動(dòng)分離并制造出一個(gè)超聲速的欠膨脹射流，射流在靠近第二錐錐面處流向下游，分離區(qū)尺寸取決于激波入射點(diǎn)位置及激波強(qiáng)度。在雙錐體的拐角上形成的分離帶和與第一個(gè)錐角上的斜激波相互作用產(chǎn)生了分離激波。

(a) Ma

注意，此風(fēng)洞實(shí)驗(yàn)將經(jīng)過Laval噴管加速的 N2作為來流氣體，由于膨脹后的N2特別稀薄所以振動(dòng)態(tài)松弛明顯，采用喉部的振動(dòng)溫度作為此來流的振動(dòng)溫度，此時(shí)氣體Tvib>>Ttr，并且在雙錐附近激波內(nèi)的流動(dòng)，依然表現(xiàn)出強(qiáng)烈的振動(dòng)態(tài)弛豫效應(yīng)。另一方面，雖然考慮了N2離解反應(yīng)，但是激波層內(nèi)生成的N原子密度極低，化學(xué)反應(yīng)過程可以忽略。

圖4展示了在CFD和DSMC方法中采用多種高溫氣體模型獲得的雙錐壁面壓強(qiáng)和熱流分布。CFD和DSMC的結(jié)果都表明，考慮振動(dòng)態(tài)激發(fā)和松弛過程條件所獲得的壁面熱流升高而壁面壓強(qiáng)無變化，因?yàn)閴簭?qiáng)主要由來流氣體動(dòng)量決定;而化學(xué)反應(yīng)十分微弱，所以化學(xué)反應(yīng)過程是否考慮對(duì)熱流和壓強(qiáng)均無影響。需要強(qiáng)調(diào)，如果考慮振動(dòng)態(tài)激發(fā)，但是認(rèn)為它總是與平動(dòng)-振動(dòng)保持平衡，會(huì)造成分離渦區(qū)域減小而二次反射波不明顯。忽略振動(dòng)態(tài)激發(fā)會(huì)造成第一錐面上的熱流預(yù)估值偏小，這是由于雙錐流動(dòng)實(shí)驗(yàn)中，來流的振動(dòng)能遠(yuǎn)高于平動(dòng)能。CFD和DSMC在不同設(shè)置下獲得的壁面壓強(qiáng)和熱流的分布，具有很好的一致性，且都與實(shí)驗(yàn)值符合。這說明，CFD和DSMC方法分別從宏觀和微觀發(fā)展了對(duì)熱化學(xué)非平衡現(xiàn)象的描述方法，均能很好地描述此復(fù)雜流動(dòng)，另一方面，由于此流動(dòng)比較稀薄，所以DSMC方法是適用的。后續(xù)的中焓雙楔-1和高焓雙楔-2流動(dòng)十分稠密，而DSMC方法要求網(wǎng)格尺寸要嚴(yán)格小于局地的分子平均自由程，受限于計(jì)算能力，后續(xù)模擬只使用了CFD方法。

(a) Pressure

2.2 中焓雙楔-1流動(dòng)

圖5展示了CFD方法中是否考慮振動(dòng)態(tài)激發(fā)與松弛以及N2的離解過程(CFD-noVib，CFD-Vibeq，CFD-Vibneq和CFD-Vibneq+React)對(duì)中焓雙楔-1溫度場(chǎng)預(yù)測(cè)的影響?？紤]激發(fā)振動(dòng)態(tài)可以明顯降低流場(chǎng)內(nèi)溫度的最大值，減小激波層厚度、分離區(qū)尺寸、三叉點(diǎn)位置等流場(chǎng)特征，并且CFD-Vibneq獲得的平動(dòng)溫度介于CFD-noVib和CFD-Vibeq的溫度場(chǎng)之間，而所獲得的振動(dòng)溫度明顯低于平動(dòng)溫度，在兩個(gè)楔面附近的振動(dòng)溫度都不高于3 000 K，而平動(dòng)溫度主要介于3 000～6 000 K之間。

圖6展示了考慮高溫條件下振動(dòng)態(tài)激發(fā)松弛過程以及N2的離解過程對(duì)雙楔-1壁面氣動(dòng)參數(shù)分布的影響?？梢钥闯鯟FD-noVib獲得的分離渦尺寸最大，CFD-Vibneq的結(jié)果是介于CFD-Vibeq和CFD-noVib之間，而CFD-Vibneq+React與CFD-Vibneq相近。仿真獲得的熱流僅在第一楔面開始部分與試驗(yàn)結(jié)果一致，之后的分離區(qū)和第二楔面上都存在明顯差別，本文的仿真結(jié)果與文獻(xiàn)[28]中Komives等的仿真結(jié)果基本相同，這種偏差很可能是由三維效應(yīng)造成的。

(a) Temperature of CFD-noVib

2.3 高焓雙楔-2流動(dòng)

圖7展示了CFD方法中是否考慮振動(dòng)態(tài)激發(fā)與松弛以及N2的離解過程(CFD-noVib，CFD-Vibeq，CFD-Vibneq和CFD-Vibneq+React)對(duì)高焓雙楔-2溫度場(chǎng)預(yù)測(cè)的影響，和前面算例一樣未使用湍流模型。當(dāng)忽略了N2振動(dòng)態(tài)激發(fā)時(shí)，流場(chǎng)溫度最大值將近20 000 K，考慮振動(dòng)態(tài)激發(fā)的流場(chǎng)溫度最大值約為17 200 K，進(jìn)一步考慮N2離解，則流場(chǎng)溫度最大值降低至14 400 K。第一楔面附近存在強(qiáng)烈的振動(dòng)松弛效應(yīng)而化學(xué)效應(yīng)不明顯，而第二楔面附近振動(dòng)態(tài)與平動(dòng)態(tài)間基本保持平衡且N2的離解強(qiáng)烈。并且考慮振動(dòng)激發(fā)和離解反應(yīng)所獲得的斜激波層更薄、三叉點(diǎn)前移且分離區(qū)減小。

圖6 雙楔-1壁面熱流分布Fig.6 Wall heat flow distribution of the double wedge-1

(a) Temperature of CFD-noVib

圖8展示了考慮高溫條件下振動(dòng)態(tài)激發(fā)松弛過程以及N2的離解過程對(duì)雙楔-2壁面氣動(dòng)參數(shù)分布的影響?？梢钥闯鯟FD-noVib獲得的分離渦尺寸最大，而CFD-Vibeq，CFD-Vibneq和CFD-Vibneq+React的分離渦尺寸較小且相近。

圖8 雙楔-2壁面熱流分布Fig.8 Wall heat flow distribution of the double wedge-2

從圖6和8可以看到，中焓雙楔-1和高焓雙楔-2流動(dòng)的CFD計(jì)算結(jié)果與實(shí)驗(yàn)均存在較明顯的偏差，分析可能存在以下原因:(1)雙楔流動(dòng)具有顯著的三維效應(yīng)，而本文直接進(jìn)行的二維模擬，這不可避免地引入誤差;(2)本文的CFD模擬均采用層流假設(shè)，高焓雙楔-2流動(dòng)的來流Re=2.9 ×105，會(huì)在第二楔面上發(fā)生轉(zhuǎn)捩過程，后續(xù)需要采用能夠預(yù)測(cè)轉(zhuǎn)捩點(diǎn)的層流-湍流模型進(jìn)一步仿真研究;(3)考慮高焓來流的穩(wěn)定性以及試驗(yàn)測(cè)量精度，實(shí)驗(yàn)結(jié)果也具有很大的不確定性。

3 結(jié)論

采用可描述熱化學(xué)非平衡過程的CFD和DSMC方法對(duì)三類雙錐/雙楔流動(dòng)進(jìn)行了數(shù)值模擬，并討論了高溫真實(shí)氣體效應(yīng)對(duì)流場(chǎng)結(jié)構(gòu)和氣動(dòng)參數(shù)的影響，得到了如下結(jié)論:

(1)隨著來流氣體焓能的增加，高溫真實(shí)氣體效應(yīng)越顯著，數(shù)值方法對(duì)熱化學(xué)非平衡過程的準(zhǔn)確模擬會(huì)影響到流場(chǎng)和氣動(dòng)分布的預(yù)測(cè)。

(2)采用量熱完全氣體模型往往會(huì)高估流場(chǎng)溫度、激波厚度，并造成拐角附近的分離渦過早分離而激波三叉點(diǎn)后移。

(3)不同工況的內(nèi)能松弛和化學(xué)反應(yīng)程度不同，比如低焓雙錐流動(dòng)中振動(dòng)能表現(xiàn)出強(qiáng)烈的松弛現(xiàn)象而化學(xué)反應(yīng)可忽略，高焓雙楔-2流動(dòng)中振動(dòng)態(tài)激發(fā)且存在明顯的化學(xué)反應(yīng)，但是振動(dòng)態(tài)松弛現(xiàn)象不明顯。在高超聲速模擬中應(yīng)注意正確應(yīng)用物理現(xiàn)象的模化模型。