考慮化學(xué)非平衡效應(yīng)的高溫湍流邊界層統(tǒng)計特性分析

劉朋欣，李 辰，孫 東，傅亞陸，袁先旭，*

(1. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室，綿陽 621000；2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000)

0 引言

Morkovin假設(shè)[1]是否成立是可壓縮湍流邊界層研究的一個重要部分。該假設(shè)指出，對于中等馬赫數(shù)的自由來流（Ma≤5），可壓縮湍流與不可壓縮湍流的區(qū)別可以通過引入流體特性的平均變化來消除。這個理論是van Direst變換的基礎(chǔ)，可壓縮邊界層內(nèi)的流向速度經(jīng)過變換后，其分布能與不可壓縮邊界層較好地符合在一起。較多的試驗和數(shù)值模擬結(jié)果[2-3]都證實了van Direst變換后邊界層速度型分布的一致性，且在高超聲速湍流邊界層中也基本成立[4-7]。另外一個在不可壓縮湍流邊界層中常用的分析是強(qiáng)雷諾比擬關(guān)系（strong Reynolds analogy，SRA），它表征了溫度脈動與速度脈動之間的關(guān)系。SRA關(guān)系式也已經(jīng)在較為寬泛的馬赫數(shù)范圍和不同的壁面溫度下得到了評估和驗證[7-11]，并發(fā)展了一些修正的SRA關(guān)系式，比如GSRA[9]、RSRA[10]、HSRA[11]。

高超聲速飛行器在中低空以極高馬赫數(shù)飛行時，頭部會產(chǎn)生較強(qiáng)的激波并壓縮來流空氣，使其溫度急劇升高，導(dǎo)致化學(xué)非平衡效應(yīng)的出現(xiàn)。 基于量熱完全氣體假設(shè)的一些湍流邊界層經(jīng)典假設(shè)和關(guān)系式在高溫化學(xué)非平衡湍流邊界層中是否仍然成立，有待于進(jìn)一步的驗證。Duan等[12]研究了焓值對高超聲速湍流邊界層的影響，在高焓時考慮空氣化學(xué)反應(yīng)模型，而低焓時則使用量熱完全氣體模型。他們發(fā)現(xiàn)，低焓流動中大部分的尺度關(guān)系，比如邊界層平均速度剖面、SRA關(guān)系式，通過消除量熱完全氣體假設(shè)后仍然能夠成立，并提出了更通用的關(guān)系式GHSRA；高焓效應(yīng)帶來的可壓縮效應(yīng)較小，對湍流結(jié)構(gòu)的影響不明顯。Kim等[13]進(jìn)一步考慮了熱非平衡效應(yīng)的影響，對比了O2或N2單組分的完全氣體模型和化學(xué)反應(yīng)模型，研究表明，熱化學(xué)非平衡效應(yīng)下，基于Morkovin假設(shè)的相關(guān)尺度關(guān)系仍然能夠成立。

但是，目前關(guān)于熱化學(xué)非平衡湍流邊界層特性的研究還比較有限。且Duan等[12]在計算時，使用高焓化學(xué)非平衡氣體模型與低焓完全氣體模型對比，沒有在同一工況下評估兩種氣體模型之間的差異；而Kim等[13]雖然在同一工況下對比了兩種氣體模型，但采用的是O2或N2單組分工質(zhì)，沒有考慮包含多種組分的復(fù)雜反應(yīng)。有必要選擇更多的來流工況，更加廣泛地研究高溫化學(xué)非平衡湍流邊界層特性。我們之前[14]研究了高溫化學(xué)非平衡湍流邊界層的標(biāo)度律，并發(fā)現(xiàn)經(jīng)典標(biāo)度律理論[15]仍然成立。本文進(jìn)一步研究高溫化學(xué)非平衡效應(yīng)對基于Morkovin假設(shè)的湍流邊界層特性關(guān)系式的影響。

1 計算模型與方法

采用含化學(xué)反應(yīng)的三維Naiver-Stokes方程進(jìn)行直接數(shù)值模擬，具體求解過程及算例驗證可參考文獻(xiàn)[14，16]。無黏通量的離散采用七階WENO-Z格式[17]，黏性通量的離散采用四階中心差分格式，時間推進(jìn)采用顯式的三階TVD（Total Variation Diminishing）型Runge-Kutta法?？諝饣瘜W(xué)反應(yīng)模型采用5組分（N2，O2，N，O，NO）6基元反應(yīng)模型[18]，只考慮空氣組分的離解。

計算狀態(tài)示意圖見圖1。在30 km高度、以馬赫數(shù)20進(jìn)行飛行的楔形體，會在頭部產(chǎn)生一道斜激波。選擇斜激波后的氣體狀態(tài)作為湍流邊界層外緣流動狀態(tài)：馬赫數(shù)Me=4.5， 溫度Te=3400K，密度ρe=0.10248kg/m3。此時溫度已經(jīng)足夠激發(fā)空氣發(fā)生化學(xué)反應(yīng)，產(chǎn)生非平衡效應(yīng)。

圖1 計算模型狀態(tài)示意圖Fig. 1 Schematic of computational condition

在考慮高溫湍流邊界層中的化學(xué)非平衡效應(yīng)時，來流組分的質(zhì)量分?jǐn)?shù)設(shè)置為f(N2) = 0.7325，f(O2) =0.2675，記為Multi算例。作為對比，本文還使用了同一來流狀態(tài)的完全氣體模型進(jìn)行計算，以評估完全氣體模型的適用性及化學(xué)非平衡效應(yīng)的影響，并記為Perf算例。

初場采用相同來流狀態(tài)下SST-RANS[19]（Shear Stress Transfer-Reynolds Averaged Naiver-Stokes）計 算結(jié)果的剖面生成。湍流入口脈動的生成采用數(shù)字濾波合成湍流方法[20]；出口為特征邊界條件；上邊界設(shè)置為初始值，并保持不變；壁面設(shè)置為黏性等溫，壁溫Tw等 于Te；展向為周期邊界。

邊界層厚度和雷諾數(shù)見表1。此時的邊界層名義厚度 δ約為5 mm，基于動量厚度的雷諾數(shù)Reθ約為2600。分析所采用流向位置如圖2中黑色線所示，且如無特殊說明，下文中分析默認(rèn)基于此流向截面處的流場。并在下文的表述中約定：qˉ代 表變量q的雷諾平均，對應(yīng)的脈動量為（ 或 ）表示其Favre平均，對應(yīng)的脈動量為。

圖2 Multi算例瞬時三維流場和展向中間截面參數(shù)分布Fig. 2 Instantaneous 3D flow field and parameter contours in the spanwise middle plane for the Multi case

表1 邊界層厚度、雷諾數(shù)和網(wǎng)格分辨率參數(shù)[14]Table 1 Parameters for the boundary layer thickness, Reynold number and grid resolution[14]

計算域設(shè)置為流向約20δ，法向約4δ，展向約2δ。湍流充分發(fā)展段的網(wǎng)格分辨率見表1。在之前的研究[14]中，進(jìn)行了同等工況下的網(wǎng)格無關(guān)性驗證，最終采用的網(wǎng)格為901×161×301（流向×法向×展向），可以滿足直接數(shù)值模擬的要求。

圖2（a）給出了瞬時流場的Q準(zhǔn)則等值面（Qcr=0.05），并用O組分質(zhì)量分?jǐn)?shù)著色。可見經(jīng)過一段距離后，湍流邊界層得到充分發(fā)展。O的生成主要集中在壁面附近溫度較高的區(qū)域。展向中間截面的瞬時溫度和O組分質(zhì)量分?jǐn)?shù)分布見圖2（b）和圖2（c），邊界層的溫度在局部已經(jīng)可以達(dá)到5000 K以上，且壁面為3400 K的冷壁，較大的溫度梯度和高溫來流條件會使得壁面產(chǎn)生較大熱流，使飛行器面臨惡劣的熱環(huán)境；在湍流充分發(fā)展區(qū)，已經(jīng)有大量的O2分子發(fā)生離解反應(yīng)并生成O原子，O原子的質(zhì)量分?jǐn)?shù)最大已經(jīng)超過0.12，這表明此時發(fā)生了較強(qiáng)的化學(xué)非平衡效應(yīng)。圖2（d）給出了展向中間截面的瞬時無量綱密度梯度大小，也可以清楚地看到，本文分析所采用的流向位置處的湍流邊界層已經(jīng)充分發(fā)展。

2 結(jié)果與討論

2.1 平均量分析

課題組前期工作中[14]分析了相同計算設(shè)置下高溫非平衡湍流邊界層的平均速度剖面、溫度剖面和雷諾應(yīng)力強(qiáng)度。這里只給出了平均溫度剖面，見圖3，以方便后續(xù)分析討論。研究發(fā)現(xiàn)平均速度剖面仍然存在明顯的對數(shù)區(qū)，對數(shù)區(qū)內(nèi)斜率基本保持不變（為1/k，k為0.41），截距有所升高C= 6.2。相較于完全氣體模型，考慮空氣化學(xué)非平衡效應(yīng)會使得邊界層內(nèi)溫度顯著降低，這是由于空氣發(fā)生了較強(qiáng)的吸熱離解反應(yīng)，使得組分內(nèi)能轉(zhuǎn)換為化學(xué)能。采用摩擦速度無量綱化的雷諾應(yīng)力分布與文獻(xiàn)中超聲速湍流邊界層[21-22]、高焓湍流邊界層[12]結(jié)果一致。這說明采用歸一化方法后，高溫化學(xué)非平衡湍流邊界層中的平均速度剖面和雷諾應(yīng)力仍然與完全氣體具有相似的分布。

圖3 平均溫度的法向剖面[14]Fig. 3 Wall-normal profiles of the averaged temperature[14]

2.2 Walz關(guān)系式

對于零壓力梯度的平板邊界層，Walz提出的修正Crocco關(guān)系式[23]（稱為Walz關(guān)系式）來建立平均溫度和平均速度間的關(guān)聯(lián)，具體表達(dá)式為：

其 中，Taw為 恢 復(fù) 溫 度恢復(fù)因子r= 0.9。

圖4給出了DNS得到的平均速度-平均溫度的關(guān)聯(lián)曲線與Walz關(guān)系式之間的對比，可以看到兩者存在很大的差異，尤其是對于空氣化學(xué)反應(yīng)模型得到的結(jié)果差異更大，這說明高溫化學(xué)非平衡效應(yīng)對平均溫度-速度的關(guān)聯(lián)性有相當(dāng)大的影響。這可能是由于高溫氣體效應(yīng)的吸熱反應(yīng)對平均溫度分布有明顯的降低作用，由此導(dǎo)致邊界層內(nèi)溫度分布發(fā)生顯著改變。

圖4 Walz關(guān)系式Fig. 4 Walz relationship

為了消除熱力學(xué)狀態(tài)、化學(xué)反應(yīng)對溫度-速度關(guān)系的明顯依賴，可以引入了一個無量綱參數(shù)—恢復(fù)焓并建立恢復(fù)焓與速度之間的關(guān)系。定義為：其中恢復(fù)因子仍設(shè)置為r= 0.9。在壁面處，在邊界層邊緣y=δ， 有則。圖5給出了與 速度之間的關(guān)系?？梢钥吹脚c的值接近，并且與文獻(xiàn)[12]符合較好。這說明存在如下關(guān)系式，或者

圖5 無量綱恢復(fù)焓與速度之間的關(guān)系Fig. 5 Relation between the non-dimensional recovery enthalpy and the velocity

經(jīng)過對DNS數(shù)據(jù)的擬合可得：

將式（3）代入式（2），并給出修改后的平均靜焓與速度之間的關(guān)系（圖6）。可見采用式（2）得到的曲線與DNS數(shù)據(jù)吻合較好。這說明在邊界層內(nèi)采用恢復(fù)焓，并使用擬合得到的關(guān)系式去修正平均靜焓與速度的關(guān)系式是有效的，可以消除自由來流的馬赫數(shù)、壁溫、化學(xué)反應(yīng)等因素的影響，且對不同工況具有一定的普適性。

圖6 修正后的～ 關(guān)系式Fig. 6 Modified ～relation

2.3 強(qiáng)雷諾比擬分析

可壓縮邊界層內(nèi)各參數(shù)脈動量之間的關(guān)系也很重要，這可以通過SRA來表征。MorKovin[1]在1962年提出了5個強(qiáng)雷諾比擬關(guān)系式，其中4個為：

對應(yīng)于上述三種改進(jìn)，c=1.0、c=1.34、c=Prt。 其中Ma為當(dāng)?shù)伛R赫數(shù)，為 總 溫。前 期 的 研 究中 顯 示[2,7,24]，HSRA對低焓情況下的絕熱或非絕熱湍流邊界層都有較好的適用性。以上三種修正是基于量熱完全氣體假設(shè)，另外還有一種修正關(guān)系式GHSRA[12]，其消除了完全氣體假設(shè)并適用于考慮化學(xué)非平衡的工況，具體為：

圖7（b）給出了基于空氣化學(xué)反應(yīng)模型算例并采用四種不同SRA修正關(guān)系式的分布，可以看到考慮了化學(xué)非平衡過程的GHSRA的值基本上在1附近，表現(xiàn)最好?？拷诿娓浇耐蛔兒头柕母淖兪怯捎谕蛔兏浇臏囟忍荻鹊扔?引起的，見平均溫度分布圖（3）。圖7（a）還對比了采用原始SRA公式（4）與GHSRA公式（9）時，兩種不同氣體模型的溫度脈動-流向速度脈動關(guān)系式的分布，可以看到兩種氣體模型下GHSRA的值都基本上在1附近，這也說明此式在消除量熱完全氣體假設(shè)時的有效性。

圖7 溫度脈動與流向速度脈動的強(qiáng)雷諾比擬關(guān)系Fig. 7 Strong Reynolds analogy relation between the temperature fluctuation and the streamwise velocity fluctuation

圖8 給出了流向速度脈動與溫度脈動之間的關(guān)系式Ru′′T′′、法向速度脈動與溫度脈動之間的關(guān)系式Rv′′T′′、流向速度脈動與法向速度脈動之間的關(guān)系式Ru′′v′′的 分布?？梢钥吹皆谶吔鐚油鈪^(qū)，u′′和T′′并不是完美的反相關(guān)性，Ru′′T′′大部分為?0.6～?0.7之間，這與文獻(xiàn)[2，12，24-26]計算結(jié)果一致。式（6）表明Ru′′v′′與Rv′′T′′呈現(xiàn)相反的相關(guān)性，圖8中兩者的分布曲線證實了這一關(guān)系式的正確性。

圖8 不同參數(shù)脈動量之間的雷諾比擬關(guān)系Fig. 8 Strong Reynolds analogy relation between the fluctuations of different parameters

式（7）給出了湍流Prandtl數(shù)的定義，它表征了平均運動中動量交換與熱交換之比。強(qiáng)雷諾比擬中認(rèn)為該值近似于1。圖9中給出了兩種氣體模型得到的Prt的分布，可以看到Prt約在0.8附近波動。這與Pirozzoli等[3]、Duan等[8]和Li等[27]計算較高來流馬赫數(shù)的工況結(jié)果一致。另外，圖9中還給出了表征湍流質(zhì)量擴(kuò)散的特征參數(shù)Prρ的分布，其定義為：

圖9 Prandtl數(shù)的法向分布Fig. 9 Wall-normal distributions of the Prandtl number

Pirozzoli等[3]認(rèn)為Prt≈Prρ，即在Morkovin假設(shè)下，湍流邊界層中動量交換、熱量交換和質(zhì)量交換過程處于近似相當(dāng)?shù)臓顟B(tài)。圖9的結(jié)果表明，在邊界層中的大部分區(qū)域，采用完全氣體模型時，Prρ的大小與Prt基本相當(dāng)。但當(dāng)考慮高溫化學(xué)非平衡效應(yīng)而采用空氣化學(xué)反應(yīng)模型時，Prρ與Prt之間的差異稍大；且越靠近邊界層外層，兩者的差異越明顯，說明三種交換過程的相對大小發(fā)生了變化。

以上分析表明，當(dāng)前計算的工況下的高溫化學(xué)非平衡湍流邊界層中，強(qiáng)雷諾比擬關(guān)系式及其修正基本上仍然能夠成立。

2.4 湍動能生成耗散機(jī)制

湍動能輸運方程常被用來分析湍動能生成、耗散等機(jī)制。湍動能的定義為：

其輸運方程為：

圖10給出了湍動能輸運方程各項沿壁面法向的分布，并采用壁面參數(shù)進(jìn)行歸一化，其中。從圖10（a）可以看到，高溫化學(xué)非平衡湍流邊界層內(nèi)湍動能的生成項P、輸運項T、擴(kuò)散項D和黏性耗散項ε占主導(dǎo)地位，其他各項的值都比較小。圖10（b）對比了兩種氣體模型的結(jié)果。在采用壁面尺度歸一化的法向坐標(biāo)y+下，考慮化學(xué)非平衡效應(yīng)時得到的生成項P、輸運項T、擴(kuò)散項D的峰值更靠近于壁面。這可能是由于兩種氣體模型計算得到的壁面附近參數(shù)不同所導(dǎo)致的。

圖10 壁面尺度下湍動能輸運方程各項的法向分布Fig. 10 Wall-normal distributions of the turbulent kinetic energy budget terms scaled in wall units

圖11 半當(dāng)?shù)爻叨认峦膭幽茌斶\方程各項的法向分布Fig. 11 Wall-normal distributions of the turbulent kinetic energy budget terms scaled in semi-local units

2.5 可壓縮效應(yīng)

本文計算的工況馬赫數(shù)較高，且考慮了化學(xué)非平衡效應(yīng)，因此有必要分析可壓縮效應(yīng)的影響。湍流馬赫數(shù)Mt常用來表征可壓縮效應(yīng)的大小[26-27]，其定義為：

且通常認(rèn)為Mt＞ 0.3時就要考慮可壓縮效應(yīng)對湍流特性的影響。圖12給出了Mt在邊界層內(nèi)的分布。可以看到，在近壁區(qū)域，兩種氣體模型得到的Mt均大于0.3，且空氣反應(yīng)模型的值略大。這可能是由于化學(xué)反應(yīng)使得Multi算例中氣體的比熱降低、聲速降低導(dǎo)致。

圖12 湍流馬赫數(shù)的法向分布Fig. 12 Wall-normal distributions of the turbulent Mach number

進(jìn)一步可以通過脹量項來研究可壓縮效應(yīng)的影響，包括壓力-脹量項 Πd和膨脹-耗散項 εd，這兩項都是由于可壓縮性導(dǎo)致速度散度不為零引起的。圖13給出了壓力做功項 Π中三個分項 Πp、 Πt和 Πd沿壁面法向的分布，并采用生成項P進(jìn)行歸一化?？梢钥吹綁毫ο嚓P(guān)項 Π主要來自于壓力輸運 Πp的作用，其余兩項所占比例很小。在邊界層中，壓力-膨脹項的大小不足湍動能生成項的3%。膨脹-耗散項與螺旋耗散項之比 εd/εs，可以用來表征黏性耗散中壓縮性部分與不可壓縮部分的比例。圖14給出了兩者之比的分布，可以看到在邊界層中 0

圖13 壓力相關(guān)項的法向分布Fig. 13 Wall-normal distributions of the pressure-related terms

圖14 膨脹-耗散與螺旋耗散之比的法向分布Fig. 14 Wall-normal distributions of the dilatational to solenoidal dissipation ratio

3 結(jié)論

本文對空間發(fā)展的高超聲速高溫湍流邊界層進(jìn)行了直接數(shù)值模擬，結(jié)果表明：

1）Walz關(guān)系式給出的平均溫度-速度關(guān)系與DNS的結(jié)果相差較大；采用無量綱恢復(fù)焓與速度之間的關(guān)系式能夠消除馬赫數(shù)、化學(xué)反應(yīng)等因素的影響，與DNS結(jié)果符合較好。

2）GHSRA可以較好地描述溫度脈動與速度脈動之間的關(guān)系，強(qiáng)雷諾比擬關(guān)系式及其修正在高溫化學(xué)非平衡湍流邊界層中基本上成立。

3）湍動能的生成項、輸運項、擴(kuò)散項和黏性耗散項在湍動能輸運過程中占主導(dǎo)地位，且使用半當(dāng)?shù)爻叨葰w一化各輸運項和壁面法向高度時，不同工況下湍動能輸運方程各項的分布能夠較好地符合。壓力-脹量項和膨脹-耗散項的值較小，化學(xué)非平衡效應(yīng)及高馬赫數(shù)效應(yīng)引起的可壓縮效應(yīng)有限。