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

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

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

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

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

對于非定常可壓縮湍流，雷諾應力采用線性渦黏性模型計算時，質量平均下的平均質量、動量和能量方程分別為

(1)

(2)

(3)

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

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

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

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

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

圖1 壓縮拐角實驗模型示意圖[15]Fig.1 Test model of compression corner experiment[15]

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

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

圖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 計算結果與分析

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

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

3.1.1 單項影響規(guī)律

圖4給出了采用no-A、no-B和no-C組合(忽略任意一個湍動能項)與采用All組合計算得到的物面壓力、熱流和摩阻結果。從圖4中可以看到，當湍動能項A被單獨忽略后，干擾區(qū)壓力峰值變化不大，但熱流峰值有一個較為明顯的增加(大約20%)；另一個更為顯著的影響則是流動分離區(qū)，在A項被忽略后，同All組合相比，其范圍縮小了近60%。當湍動能項B被單獨忽略時，壓力和熱流峰值的變化均較小，但是流動分離區(qū)也同樣發(fā)生了非常明顯的變化；忽略了平均總能量定義中的B項會大大加速流動分離的發(fā)生(與A項被忽略的情形剛好相反)，此時整個分離區(qū)比All組合大了近30%。而當湍動能項C被單獨忽略時，與B項被單獨忽略時變化情況相似：物面壓力峰值和熱流峰值都幾乎不受影響，但分離區(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給出了忽略任意兩個湍動能相關項后(no-A_B、no-A_C和no-B_C組合)計算得到的物面壓力、熱流和摩阻與包含全部3個湍動能項時的對比情況?？梢钥吹?，4種情況得到的壓力峰值基本相同，而干擾區(qū)熱流峰值在A項被忽略的條件下，不論B項或C項哪一個被同時忽略，熱流峰值都會相較于完整包含3個湍動能項時有所升高(15%左右)，這與A項被單獨忽略時規(guī)律相似；在B項和C項被同時忽略時，熱流峰值會有所下降(10%左右)。從物面摩阻系數(shù)的分布來看，流動分離區(qū)受A項的影響最大：當這一項被忽略后(不管其他兩項哪一個被一同忽略)，分離區(qū)范圍都會縮小，只是在包含B項時(no-A_C組合)分離區(qū)減小的程度相對弱一些。當B項和C項被同時忽略時，分離區(qū)的變化情況與B項或C項被單獨忽略時的規(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個湍動能相關項都被忽略的方程(no-A_B_C組合)與All組合計算所得到的物面壓力、熱流和摩阻的對比情況。不難發(fā)現(xiàn)，當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個)對物面壓力峰值的影響很小，相應的計算結果與包含全部3個湍動能項時基本一致。干擾區(qū)熱流峰值受湍動能的影響相較于壓力峰值的變化要更明顯得多，尤其是當湍動能項A被忽略時，熱流峰值會上升12%～20%。受湍動能項影響最明顯是流動分離，表4列出了8種不同湍動能組合條件下，分離點、再附點和分離區(qū)范圍。從表4中可以看出，只要湍動能項A被忽略，同包含全部3個湍動能項相比，流動分離都會滯后發(fā)生；這種影響在A項被單獨忽略時最為明顯，整個分離區(qū)的范圍可以比不省略時小近60%。當包含A項時(完整的平均動量方程)，忽略湍動能項B或C都會加速流動分離的發(fā)生，但是C項單獨作用的影響相對較小，分離區(qū)范圍與不省略時增加8%左右。同時忽略B項和C項，對流動分離的加速效果最為顯著，其分離區(qū)范圍可以達到不忽略時(包含全部3個湍動能項)的近1.4倍。

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

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

3.2 湍動能項作用機理

3.2.1 湍動能項A作用機理

(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作用機理

(5)

(6)

(7)

(8)

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

圖10 All、no-A和no-B湍動能項組合在分離點附近有效壓力場流向壓力梯度分布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作用機理

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

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

4 結 論

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

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

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

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

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

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