張昊元，孫東，邱波，朱言旦,王安齡,*

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

以高超聲速在大氣層內(nèi)飛行是未來飛行器發(fā)展的一個重要方向，相關(guān)技術(shù)吸引了眾多研究機(jī)構(gòu)和學(xué)者的持續(xù)關(guān)注。激波/邊界層干擾現(xiàn)象廣泛存在于這一類飛行器的流場中，其通常發(fā)生在飛行器的不同組件之間，常常引起氣動力/熱等關(guān)鍵參數(shù)的顯著改變。對激波/邊界層干擾流動的準(zhǔn)確預(yù)測，是整個飛行器氣動設(shè)計中的重要一環(huán)。

為此，本文針對目前對湍動能在RANS方程中，特別是在涉及激波/邊界層干擾的模擬中作用不清的現(xiàn)狀，采用基于OpenFOAM平臺開發(fā)的求解器，針對=7.05的柱-裙實驗?zāi)Ｐ烷_展了數(shù)值模擬計算，對湍動能的影響規(guī)律進(jìn)行了研究，并系統(tǒng)分析了造成這些影響的物理機(jī)制。

1 控制方程與湍動能項組合方式

對于非定?？蓧嚎s湍流，雷諾應(yīng)力采用線性渦黏性模型計算時，質(zhì)量平均下的平均質(zhì)量、動量和能量方程分別為

(1)

(2)

(3)

從式(2)和式(3)可以發(fā)現(xiàn)，有3類與湍動能相關(guān)的項。其中A項是渦黏性模型雷諾應(yīng)力定義式中的一部分，B項是平均總能量定義中的湍動能分量，C項是分子擴(kuò)散與湍流輸運(yùn)項(采用簡單梯度方法?；?。當(dāng)A～C項都被省略后，就是引言中提到的CFL3D或rhoCentralFoam中處理RANS模擬時求解的方程平均量輸運(yùn)方程。這里值得指出的是，在CFL3D中采用非線性渦黏性模型進(jìn)行計算時，求解的RANS方程中與采用線性渦黏性模型時是有所不同的，其中只有湍動能項B被忽略了，而湍動能項A和C都保留了。

為了系統(tǒng)研究湍動能項A、B和C對激波/邊界層干擾流動的影響規(guī)律，提出了8種不同的湍動能相關(guān)項組合方式(其中有3種與文獻(xiàn)[8]中相同)。具體組合信息由表1給出，其中CFL3D軟件在處理線性模型時采用的就是“no-A_B_C”這一組合方式,而處理非線性模型時則變成了“no-B”這一種組合。由此可見，當(dāng)上述湍動能項(一個或其中幾個)對預(yù)測結(jié)果有顯著影響時，采用含有不同湍動能相關(guān)項組合的RANS方程對不同湍流模型的性能與表現(xiàn)進(jìn)行評價時，是缺乏相同的比較基礎(chǔ)的(求解了不同形式的平均流場輸運(yùn)方程)，因此采用這種方式得出的結(jié)論也就包含了很大的不確定性。

表1 湍動能項組合方式Table 1 Combinations of turbulent kinetic energy terms

2 數(shù)值方法與算例設(shè)置

采用有限體積方法在三維多邊形網(wǎng)格上對文中涉及到的所有控制方程進(jìn)行離散求解，數(shù)值計算均在OpenFOAM開源平臺上實現(xiàn)?；贠penFOAM內(nèi)置的rhoCentralFoam求解器開發(fā)了新的基于密度的求解器，并包含了表1中列出的全部湍動能項組合。數(shù)值計算中所有對流項(包括平均量和湍流變量的輸運(yùn)方程)都采用AUSMPW+格式離散，在交界面上原始變量的重構(gòu)采用帶有Van-Albada限制器的二階TVD方法實現(xiàn)，而擴(kuò)散項計算均采用中心差分格式，時間推進(jìn)采用基于局部時間步長的一階隱式歐拉方法，為了保證數(shù)值穩(wěn)定性，所有計算的CFL數(shù)控制在5以內(nèi)。湍流模型采用Menter提出的-SST模型，其中湍流變量在物面和來流的邊界條件均按照文獻(xiàn)[14]中的推薦值設(shè)置。

圖1 壓縮拐角實驗?zāi)Ｐ褪疽鈭D[15]Fig.1 Test model of compression corner experiment[15]

表2 壓縮拐角算例無干擾邊界層參數(shù)

為檢驗計算網(wǎng)格對數(shù)值計算結(jié)果的影響，我們采用了4套網(wǎng)格(具體信息如表3所示)對網(wǎng)格收斂性進(jìn)行了分析。圖3給出了使用4套網(wǎng)格采用包含全部3個湍動能項(All組合)的RANS方程計算得到的物面壓力和摩阻結(jié)果(圖中為沿流向的坐標(biāo)，“0”點(diǎn)位于壓縮拐角處)。從中可以看出，隨著網(wǎng)格加密，流動分離區(qū)范圍會逐漸增大直至收斂，其中由Grid_fine(300×300)給出的結(jié)果與Grid_extra fine(400×400)給出的結(jié)果幾乎完全相同，說明數(shù)值計算達(dá)到了很高的網(wǎng)格收斂水平。因此，后續(xù)所有數(shù)值計算均采用Grid_Fine網(wǎng)格進(jìn)行。

圖2 壓縮拐角算例網(wǎng)格示意圖Fig.2 Grid for compression corner case

表3 網(wǎng)格信息Table 3 Grid information

圖3 不同網(wǎng)格物面壓力和摩阻對比Fig.3 Comparison of predicted wall pressure and skin friction by different grids

3 計算結(jié)果與分析

3.1 湍動能項的影響規(guī)律

這里將8種不同的湍動能項組合分成3類情況進(jìn)行比較：任意一項被單獨(dú)忽略、任意兩項被同時忽略以及三項同時被忽略。分析的重點(diǎn)是對湍動能項被忽略后，物面參數(shù)與流場結(jié)構(gòu)相較于不忽略時的變化特征。

3.1.1 單項影響規(guī)律

圖4給出了采用no-A、no-B和no-C組合(忽略任意一個湍動能項)與采用All組合計算得到的物面壓力、熱流和摩阻結(jié)果。從圖4中可以看到，當(dāng)湍動能項A被單獨(dú)忽略后，干擾區(qū)壓力峰值變化不大，但熱流峰值有一個較為明顯的增加(大約20%)；另一個更為顯著的影響則是流動分離區(qū)，在A項被忽略后，同All組合相比，其范圍縮小了近60%。當(dāng)湍動能項B被單獨(dú)忽略時，壓力和熱流峰值的變化均較小，但是流動分離區(qū)也同樣發(fā)生了非常明顯的變化；忽略了平均總能量定義中的B項會大大加速流動分離的發(fā)生(與A項被忽略的情形剛好相反)，此時整個分離區(qū)比All組合大了近30%。而當(dāng)湍動能項C被單獨(dú)忽略時，與B項被單獨(dú)忽略時變化情況相似：物面壓力峰值和熱流峰值都幾乎不受影響，但分離區(qū)范圍有所增大。不過從分離區(qū)范圍的變化情況可知，C項對流動分離帶來的影響要明顯小于B項。

圖4 忽略單個湍動能項與包含全部湍動能項計算得到的物面壓力、熱流和摩阻對比Fig.4 Comparison of predicted wall pressure, heat flux and skin friction between case with a single TKE term neglected and that with all TKE terms included

3.1.2 雙項影響規(guī)律

圖5給出了忽略任意兩個湍動能相關(guān)項后(no-A_B、no-A_C和no-B_C組合)計算得到的物面壓力、熱流和摩阻與包含全部3個湍動能項時的對比情況?？梢钥吹剑?種情況得到的壓力峰值基本相同，而干擾區(qū)熱流峰值在A項被忽略的條件下，不論B項或C項哪一個被同時忽略，熱流峰值都會相較于完整包含3個湍動能項時有所升高(15%左右)，這與A項被單獨(dú)忽略時規(guī)律相似；在B項和C項被同時忽略時，熱流峰值會有所下降(10%左右)。從物面摩阻系數(shù)的分布來看，流動分離區(qū)受A項的影響最大：當(dāng)這一項被忽略后(不管其他兩項哪一個被一同忽略)，分離區(qū)范圍都會縮小，只是在包含B項時(no-A_C組合)分離區(qū)減小的程度相對弱一些。當(dāng)B項和C項被同時忽略時，分離區(qū)的變化情況與B項或C項被單獨(dú)忽略時的規(guī)律相同，不過分離區(qū)范圍產(chǎn)生了更大幅度的增大。

圖5 忽略兩個湍動能項與包含全部湍動能項計算得到的物面壓力、熱流和摩阻對比Fig.5 Comparison of predicted wall pressure, heat flux and skin friction between case with double TKE terms neglected and that with all TKE terms included

3.1.3 三項影響規(guī)律

圖6給出了采用3個湍動能相關(guān)項都被忽略的方程(no-A_B_C組合)與All組合計算所得到的物面壓力、熱流和摩阻的對比情況。不難發(fā)現(xiàn)，當(dāng)3個湍動能都被忽略時，干擾區(qū)壓力峰值基本不受影響，這說明了這3個湍動能項對物面壓力峰值的影響非常小。而干擾區(qū)熱流峰值有所上升(12%左右)，分離區(qū)范圍則與包含3個湍動能項時相比依舊呈顯著減小的趨勢。對比B項和C項同時被忽略時(no-B_C組合)流動分離區(qū)的變化情況(顯著增大)，可以發(fā)現(xiàn)A項對流動分離的影響比B項和C項同時作用的影響還要大。

圖6 忽略3個湍動能項與包含全部湍動能項計算得到的物面壓力、熱流和摩阻對比Fig.6 Comparison of predicted wall pressure, heat flux and skin friction between case with triple TKE terms neglected and that with all TKE terms included

從上述物面參數(shù)(壓力、熱流和摩阻)的變化規(guī)律可以看出，忽略部分或者全部的湍動能項(本文中列出的3個)對物面壓力峰值的影響很小，相應(yīng)的計算結(jié)果與包含全部3個湍動能項時基本一致。干擾區(qū)熱流峰值受湍動能的影響相較于壓力峰值的變化要更明顯得多，尤其是當(dāng)湍動能項A被忽略時，熱流峰值會上升12%～20%。受湍動能項影響最明顯是流動分離，表4列出了8種不同湍動能組合條件下，分離點(diǎn)、再附點(diǎn)和分離區(qū)范圍。從表4中可以看出，只要湍動能項A被忽略，同包含全部3個湍動能項相比，流動分離都會滯后發(fā)生；這種影響在A項被單獨(dú)忽略時最為明顯，整個分離區(qū)的范圍可以比不省略時小近60%。當(dāng)包含A項時(完整的平均動量方程)，忽略湍動能項B或C都會加速流動分離的發(fā)生，但是C項單獨(dú)作用的影響相對較小，分離區(qū)范圍與不省略時增加8%左右。同時忽略B項和C項，對流動分離的加速效果最為顯著，其分離區(qū)范圍可以達(dá)到不忽略時(包含全部3個湍動能項)的近1.4倍。

表4 不同湍動能組合分離區(qū)長度對比

鑒于湍動能項對流動預(yù)測結(jié)果十分顯著的影響，分析并掌握湍動能相關(guān)項的作用機(jī)理，就十分必要了。但是在此之前，值得指出的是，從上述8種組合的計算結(jié)果與實驗結(jié)果的對比情況來看，SST模型可以較為準(zhǔn)確地給出物面壓力峰值的預(yù)測結(jié)果，但是給出的熱流峰值預(yù)測結(jié)果則嚴(yán)重高于實驗值，這與其他學(xué)者采用這一模型得到的規(guī)律相同。與此同時，對于分離區(qū)的預(yù)測，僅從這一個算例來看，采用no-A_B和no-A_B_C組合計算結(jié)果反而與測得的分離點(diǎn)位置符合得最好，而采用包含了全部湍動能項(All組合)時，分離區(qū)明顯大于實驗值。從已有的一些研究中可以發(fā)現(xiàn)，SST模型在二維/柱裙壓縮拐角類激波/邊界層干擾流動中，(由于Bradshaw數(shù)的原因)通常預(yù)測的分離區(qū)都明顯大于實驗值，這與采用All組合計算時規(guī)律是相同的。由此便可以知道，采用一個不完整的RANS方程進(jìn)行數(shù)值計算時，可以出現(xiàn)與湍流模型本身的缺點(diǎn)不足相互影響，甚至抵消缺陷的結(jié)果(當(dāng)然也可能放大這種缺陷)。這也再一次說明了采用相同的平均流動方程進(jìn)行RANS模擬計算的重要性與實際意義。

3.2 湍動能項作用機(jī)理

3.2.1 湍動能項A作用機(jī)理

(4)

圖7 包含全部湍動能項時比值分布云圖Fig.7 Contours of ratio of with all TKE terms included

圖8 采用不同湍動能項組合計算得到的平均壓力分布與DNS數(shù)據(jù)對比Fig.8 Comparison of mean pressure predicted by different TKE combinations with DNS data

圖9 no-A組合真實平均壓力場等值線分布Fig.9 Iso-lines of real mean pressure predicted by no-A combination

3.2.2 湍動能項B作用機(jī)理

(5)

(6)

(7)

(8)

以此為基礎(chǔ)，也就可以更好地解釋了在同時忽略湍動能項A和B之后(no-A_B組合)分離區(qū)依舊縮小的機(jī)制,此時的有效壓力可以近似由式(9) 計算得到：

圖10 All、no-A和no-B湍動能項組合在分離點(diǎn)附近有效壓力場流向壓力梯度分布Fig.10 Pressure gradients in streamwise direction predicted by All, no-A, and no-B TKE combinations near separation point

(9)

3.2.3 湍動能項C作用機(jī)理

對于湍動能項B和C同時被忽略的情況(no-B_C組合)，兩種相同作用的因素疊加，會使得分離區(qū)增大的程度超過單個因素起作用時的效果，這也與數(shù)值計算結(jié)果相符：no-B_C組合是所有不同組合中促使分離區(qū)增大最為顯著的一個。而湍動能組合no-A_C、no_A_B_C，則主要由于A項的影響大于其他兩個湍動能項帶來的影響，因而計算結(jié)果都呈現(xiàn)出了與單獨(dú)忽略湍動能項A時相同的變化趨勢。

從現(xiàn)有采用SST模型計算得到的物面壓力和熱流分布與實驗值的對比情況來看，數(shù)值計算的結(jié)果與實驗值相比還有較大的差距(流動分離區(qū)范圍預(yù)測偏差較大且熱流峰值被嚴(yán)重高估了)。在激波/邊界層干擾區(qū)，湍流呈現(xiàn)出顯著的各項異性特征，而提升對這種各向異性的模擬精度是改善流動預(yù)測結(jié)果的重要途徑之一。為此，在RANS模擬中可以采用更高級的湍流模型(非線性湍流模型)來實現(xiàn)對湍流各向異性特征的更為準(zhǔn)確地模擬。當(dāng)采用如非線性渦黏性模型時，在雷諾應(yīng)力的計算中包含湍動能項A對于準(zhǔn)確的反應(yīng)流場中湍流的各向異性程度和各個正應(yīng)力分量的影響都至關(guān)重要(CFL3D軟件也是采用這種方式處理非線性渦黏性模型的)。這也就意味著，采用包含全部湍動能項的RANS方程，是通過引入高級湍流模型，以期改善流動分離和物面參數(shù)預(yù)測結(jié)果的重要基礎(chǔ)。

4 結(jié) 論

本文通過對典型高超聲速激波/邊界層干擾流動的數(shù)值模擬計算，獲得了3個不同湍動能項在8種組合條件下，對流場預(yù)測結(jié)果的影響規(guī)律，并通過對平均變量輸運(yùn)方程的定性和定量分析，研究了不同湍動能項在流動預(yù)測中的作用機(jī)理，得到了以下結(jié)論：

1) 在RANS方程的中忽略雷諾應(yīng)力定義式(渦黏性模型)中湍動能項A會導(dǎo)致干擾區(qū)物面熱流峰值增大12%～20%；同時顯著推遲流動分離的發(fā)生，其分離區(qū)范圍會縮小至包含全部湍動能項時的40%左右。

2) 忽略RANS方程中平均總能量定義式中的湍動能項B和分子擴(kuò)散與湍流輸運(yùn)項C,對干擾區(qū)物面壓力和熱流影響較小，但是會加速流動分離的發(fā)生；其中B項帶來的影響要比C項嚴(yán)重得多，分離區(qū)會因B項被忽略而增大超過26%，而C項帶來的相似影響為8%左右。

4) 湍動能項C主要起著增強(qiáng)流體中湍動能擴(kuò)散的作用，因而忽略這一項會引起近壁面區(qū)域湍流黏性系數(shù)的下降，進(jìn)而產(chǎn)生加速分離的效果。3個湍動能項中，A項對流動分離的影響最強(qiáng)，B項次之，而C項最弱。

綜合上述結(jié)論可知，湍動能項在高超聲速激波/邊界層干擾流動的RANS模擬中影響十分顯著，其影響特性主要表現(xiàn)對流場中有效平均壓力場的顯著改變。因此，在涉及到受平均壓力場影響較大的流動中，例如對強(qiáng)逆壓梯度帶來的流動分離的模擬或激波誘導(dǎo)分離的臨界分離條件的模擬分析中(不限于高超聲速)，都應(yīng)該考慮采用包含全部湍動能項的RANS方程進(jìn)行模擬分析。而對于非平均壓力場主導(dǎo)流動的模擬中，如傳統(tǒng)的零壓力梯度的平板/圓柱湍流邊界層等，不同形式的RANS方程對預(yù)測結(jié)果的影響將會非常有效，這也與筆者前期的研究結(jié)論一致。

與此同時，盡管采用包含全部湍動能形式的RANS不一定總會帶來計算結(jié)果的改善(預(yù)測結(jié)果的精度還受到所采用湍流模型的影響)，但是由于忽略不同的湍動能項在不同類型流動的模擬中會產(chǎn)生不同的影響，因而采用包含全部湍動能項的RANS方程對于在相同的前提下，研究如何改進(jìn)湍流模型提升流動預(yù)測精度,具有十分重要的實際意義。