高超聲速氣動熱數值模擬的網格模式相關性研究

呂水燕，張傳俠，葉 坤，徐 健

(1.中國兵器工業試驗測試研究院， 陜西 華陰 714200； 2.西北工業大學航空學院， 西安 710072)

隨著高超聲速飛行器的進一步發展，嚴重的氣動熱是外形設計中首要解決的問題之一。關于氣動熱問題的研究方法主要有氣動熱的地面試驗研究、數值計算研究和工程估算研究[1]。相對于地面風洞試驗和實飛試驗，數值計算運行成本低、設計周期短，可以解決風洞試驗無法避免的壁面干擾和氣體物性問題。數值計算主要是求解N-S方程各種近似形式即可得到熱流分布，這種方法較為精確，在CFD計算熱流的應用中最為廣泛。美國NASA Langley研究中心開發了許多該類相關程序，最為成熟的有LAURA(Langlay Aerothermo-dynamic Upwind Relaxation Algorithm Code)、GASP(General Aerodynamic Simulation Program)、DPLR(Data-Parallel Line Relaxation Method for the Navier-Stokes Equations)等，這些軟件涵蓋了稀薄氣體效應和高溫氣體效應、湍流等[2]。美國的這些軟件已經發展應用得相當成熟，且目前氣動熱的數值計算已經開始向非結構網格技術發展[3-4]。我國氣動熱的研究起步較晚，數值計算主要通過求解NS方程預測高超聲速氣動加熱問題，國內研究以NND格式系列化研究為代表，發展了一系列提高計算精度的方法[5]。關于氣動熱的數值計算，由于對近壁的第一層高度要求很高，因此，在采用的網格模式上，一般采用結構化網格，這就在一定程度上使得外形復雜或波系干擾強烈的問題很難劃分高精度的結構網格進行數值計算求解[6-7]。閆超、潘沙等[8-9]在氣動熱的網格相關性上進行了相應的研究，而對非結構網格在氣動熱的分析上， 國內目前少有研究。

本文采用氣動熱的經典雙橢球模型為研究對象，參考文獻[10]中的試驗條件，對研究對象在8.04Ma和10.02Ma的工況下進行氣動熱數值模擬研究。計算分別采用結構網格和非結構網格，對模型表面的Y+值進行評估，計算獲得了雙橢球體表面的氣動壓力分布和熱流分布，并將兩種網格模式下的模型中心子午線上的熱流數據進行了對比分析。

1 計算模型

本文計算了文獻[10]中所述的雙橢球體標模的氣動熱特性，并與文獻中的試驗數據進行了對比分析。試驗是在Ma=8.04，Re=1.13×107/m和Ma=10.02，Re=2.20×106m兩種流場條件下進行的。本文計算針對0°攻角的計算結果進行對比分析，雙橢球模型如圖1所示。

圖1 雙橢球模型幾何圖

雙橢球體給出幾何公式如下，以保證高精度模型的生成。

試驗按照以下兩種狀態進行，如表1所示。

表1 試驗狀態

2 邊界條件推算

根據表1中的試驗狀態，推算試驗仿真邊界條件，包括：靜壓PS、靜溫TS和來流速度V。邊界條件推算采用如下公式[11]：

總溫計算公式：

馬赫數計算公式：

薩特蘭公式計算動力黏度：

薩特蘭公式計算密度：

雷諾數計算公式：

利用以上公式推導出了計算所需的參數，獲得邊界條件匯總如表2所示。

表2 根據試驗狀態推導仿真邊界

3 數值方法

結構網格劃分采用ICEM 進行，計算采用兩套網格進行計算，第1套網格單元數量為375萬，近壁層網格高度為10-5m，層高按照1∶1.1的比例向外發散，第2套網格近壁網格高度值調整為10-6m。結構網格模型如圖2(a)所示。

非結構網格同樣采用ICEM進行劃分兩套網格，近壁采用棱柱層網格劃分附面層，非結構網格的近壁網格高度同樣劃分為10-5m和10-6m，1∶1.1等比例向外發散，網格單元數554萬。非結構網格模型如圖2(b)所示。

圖2 計算模型網格

針對復雜結構模型，采用結構化網格難度大，網格生成周期也很長，難以在工程上進行應用，這就需要研究結構網格和非結構網格計算結果的差異。本文同時采用非結構網格對雙橢球體標模進行計算，與結構網格的計算結果進行對比分析，以期獲得結構網格和非結構網格計算結果差異的規律。

SST湍流模型通過混合函數將標準k-ω模型和k-ε模型結合到一起，核心思想是在近壁面利用k-ω模型的魯棒性，以捕捉到黏性底層的流動，而在主流區域則利用k-ε模型又可以避免k-ω模型對于入口湍動參數過于敏感的劣勢，SST模型在工程上應用很廣泛，本文計算采用SST湍流模型進行計算。

CFD計算格式的上風格式主要分為兩大類：以Van Leer格式為代表的FVS(Flux Vector Splitting)格式和以Roe為代表的FDS(Flux Difference Splitting)格式。針對這兩種格式各自的優缺點，又發展出AUSM[12-13](Advection Upstream Splitting Method)的混合格式。這種格式構造簡單，無矩陣運算，激波分辨率高而且穩定性好，既具有FDS格式在邊界層中解的精確性，也具有FVS格式在捕捉強間斷時的健壯性，具有非常優秀的流動計算性能。經過10多年的發展，在高超聲速的各種流動及氣動熱的數值計算研究中，獲得了廣泛的發展應用。一些商業軟件如Fluent、Overflow等在高超聲速的計算中也添加了該類格式。本文計算采用AUSM的數值離散格式。

4 不同網格模式對氣動熱的影響分析

4.1 Y+分析

熱流的仿真計算中，近壁網格密度對熱流計算結果影響很大，這已經成為共識。目前大多數討論熱流計算網格的文獻把近壁面的Y+或者網格雷諾數作為衡量仿真精度的一個非常重要的無量綱值[8-9]。該值越低，計算精度越高，但是網格尺寸也越小，網格規模越大。本文參考Fluent內核代碼中關于Y+的定義，為網格高度的特征尺寸與黏性尺度的比值，其公式為

其中:y是表征近壁網格高度的特征尺寸，對于結構網格來說y是第一層網格高度的一半，對于非結構網格來說，y第一層網格高度的1/3；μ是空氣黏性系數；ρ是空氣密度，τw壁面剪切應力。

本文從Y+入手，圖3給出了本計算中雙橢球體上下面中心子午線上的Y+的數值曲線。圖3(a)為8.04Ma工況網格子午線上的Y+值，圖3(b)為10.02Ma工況網格子午線上Y+值。由圖可以看到，計算獲得Y+數值均已控制在10以下，且結構網格和非結構網格在同等近壁網格尺度下獲得的Y+數值基本一致，在第二橢球位置之后，非結構網格的Y+數值優于結構網格。10.02Ma工況Y+數值均達到了4以下，趨勢規律與8.04Ma一致。

圖3 雙橢球體中心子午線的網格Y+值

4.2 計算結果分析

圖4以Ma=8.04攻角為0°狀態下為例，給出了流場靜壓試驗與仿真的對比。從圖4可以看到，流場的激波與試驗紋影圖吻合很好。

采用非結構網格，由于邊界層外網格的稀疏性，導致壓力云圖在激波邊界處的捕捉清晰度遠遠不如結構網格，壓力云圖的均勻性較差，流場激波區域顯得模糊。但是，流場的高壓區域的位置和結構網格的計算結果還是一致的。依照文獻[4-5]的氣動力的計算結果看，采用非結構網格，在氣動靜壓方面仿真與試驗數值吻合較好。而在氣動熱方面的差異性是本論文的重點。

壁面溫度給定280 K，采用冷壁熱流方法計算雙橢球體中心子午線上的熱流結果如圖5所示。從計算結果來看，無論是結構網格還是非結構網格，仿真計算結果的熱流分布趨勢與實驗數據均有較好的吻合度。從圖5(a)、圖5(c)可以看到，熱流順來流方向先是大幅度下降，然后趨于平緩，氣流流經第2個小橢球體時形成第2道激波，熱流在兩橢球相貫線前明顯有一個下降過程，兩種仿真邊界的結果同試驗結果一致，這表明該位置處邊界層發生了分離。小橢球處的2次激波的波后熱流迅速升高，形成了第2個熱流峰值，熱流達到峰值后逐漸下降。下表面中心線上的熱流分布如圖5(b)、圖5(d)所示，熱流經過前緣高熱流區域后逐漸下降，至模型中部趨于平緩。

圖4 流場結構試驗照片對比

圖5 雙橢球體上、下表面中心子午線熱流

由于8.04Ma流場的壓力和密度比10.02Ma的流場高出許多，參見表2，所以，8.04Ma的上表面熱流較10.02Ma高出一些。

在8.04Ma上、下表面子午線熱流曲線上，在前緣駐點位置，非結構網格的熱流數值達到了1 000 kW/m2以上，與試驗數據484 kW/m2偏差超過100%，而近壁層高0.01 mm的結構網格駐點位置的熱流計算結果為557 kW/m2與試驗結果的偏差縮小很多。這一點在10.02Ma的計算中同樣得到印證。說明非結構網格在駐點位置的熱流計算失真，這是由于在同樣計算資源和網格規模的情況下，非結構網格在前緣駐點位置也就是在熱流梯度很大，熱流變化劇烈的位置，網格密度遠遠達不到結構網格的精細度所導致的。

在離開駐點位置之后至第二橢球邊界層分離之前，在雙橢球體的上子午線試驗數據和結構網格與非結構網格的偏差量一致，兩種網格在該位置有同等計算精度。在兩橢球相貫位置，小橢球后的二次激波使得該處熱流再次激增，形成第二個熱流峰值，非結構網格在此處的模擬與試驗數據更加接近，優于結構網格的模擬精度，這是由于該處非結構獲得的Y+數值優于結構網格。隨著熱流緩慢下降，非結構網格比結構網格計算的熱流趨勢線與試驗值吻合程度更好。在下子午線，試驗數據與非結構網格偏差更小，與結構網格偏差稍有增大。

在除駐點熱流梯度很大的位置外，計算模型表面大范圍內非結構網格計算結果均優于結構網格，究其原因為本算例中，非結構網格的劃分獲得了比結構網格更優的Y+數值。這說明，在非結構網格近壁獲得更佳Y+數值的情況下，可以獲得比結構網格更好的模擬精度。

5 結論

1) 采用非結構網格計算對激波邊界的捕捉不清晰，流場激波邊界較為模糊。這是由于結構網格在外場主流區域的網格比非結構網格分布更為合理細密。

2) 本文計算熱流Y+結構網格控制在10以下，非結構網格控制在3以下。非結構網格在獲得比結構網格Y+更優數值的前提下，非結構網格在熱流曲線上比結構網格與試驗值有更高的吻合度，表面大范圍內可以獲得比結構網格更優的計算結果。

3) 非結構網格在駐點處熱流計算，誤差超過100%，這是因為在駐點處熱流梯度很大，非結構網格在駐點處的分布很難做到像結構網格一樣細密，在同等網格規模的情況下，非結構網格對高熱流梯度的位置難以準確模擬。