可壓縮流動脈動壓力數值模擬求解器HFS 研究

（中國工程物理研究院總體工程研究所，四川 綿陽 621999）

高速飛行器在飛行過程中將可能經歷高溫真實氣體效應、稀薄氣體效應、黏性干擾、邊界層轉捩和分離以及熱輻射等復雜物理化學現象[1-4]。由于實驗研究難度大，成本高，數值模擬成了重要的研究途徑之一。高速飛行中會出現各種特殊流動現象，如激波與邊界層的相互干擾、邊界層傳熱傳質、化學反應以及燒蝕。高速飛行器的內外流動以湍流及轉捩為主要特征。層流轉捩到湍流后，摩阻及熱流增加為原來的數倍。因而，無論氣動力還是熱防護設計，都需要對轉捩與湍流進行精細預測。此外，燃料混合、氣動光學以及氣動聲學都與湍流及轉捩密切相關，這就對湍流及轉捩的機理、模型及控制提出了嚴格的要求。目前數值模擬湍流的主要手段[5-8]有直接數值模擬（DNS）、大渦模擬（LES）、雷諾平均Navier-Stokes方程（RANS）以及RANS/LES 混合模型（如DES）。

直接數值模擬是目前最為準確可靠的湍流計算手段[5,9-10]，它不引入任何湍流模型，而是直接求解Navier-Stokes 方程，以獲得流場的所有時間尺度和空間尺度的細節。要捕捉全部時空演化細節，要求計算的網格尺度比流場結構的空間尺度小，同時時間步也要足夠小以捕捉相應的時間尺度，因此DNS 的網格量十分巨大。與此同時，湍流中的小尺度結構很容易被計算格式的數值耗散耗散掉，從而導致流場細節的丟失。因此，對湍流進行直接數值模擬時，要求數值格式的數值耗散要小。中心型數值格式往往具有低耗散或者無耗散，例如中心差分格式和中心型緊致格式，但是這些數值格式在計算含數值間斷問題（如激波）時，往往會在間斷附近產生數值震蕩而使計算失敗。要捕捉間斷，必須往數值格式里引入一定的數值耗散以抑制間斷附近的數值震蕩，例如使用激波捕捉格式（如TVD 格式[11]、ENO 格式[12]、WENO 格式[13]等）。高分辨率所需的低耗散與保持計算穩定的高耗散這一對矛盾很難平衡，尤其是對于高速湍流[11]。在計算流體力學（CFD）的高精度數值算法當中，基于非結構網格的高精度方法目前尚未成熟，在模擬含間斷的流動時魯棒性較差，同時計算量和內存占用較大，在實際工程中運用較少?；诮Y構網格的高精度有限差分方法具有成熟的數值方法，計算代價小，是工程中最實用的高精度方法。然而對于復雜幾何外形的流場，由于難以生成高精度的網格，實際計算中會引入網格幾何誤差，處理不慎會影響計算結果的準確性。鄧小剛等[14]指出當數值格式滿足幾何守恒律時，可以消除網格不光滑帶來的幾何誤差，其中加權緊致非線性格式（WCNS）[15-19]就是其中最典型的一種。文中求解器主要采用WCNS 格式作為空間離散格式，以確保計算高速流動中的穩定性。

綜上所述，通過數值模擬研究可壓縮流動尚存在許多技術難點。文中旨在發展一套能夠模擬復雜幾何構型的可壓縮流動的高精度求解器，并利用該求解器研究高速下鈍頭體的繞流流動。

1 數值方法

為了模擬復雜幾何構型的流場，一般采用曲線坐標系下的可壓縮Navier-Stokes 方程，其形式為：

式中對流通量分別為：

黏性通量分別為：

式中：τij是黏性應力張量；βi=μjτij+qi，qi為熱通量，黏性系數μ由Sutherland 公式確定。

為了準確捕捉激波，對流項采用顯式加權緊致非線性格式（WCNS-E）進行離散：

WCNS 格式的核心思想是，將一個大的點模板拆分為若干小的模板的線性組合，如式（5）所示。

引入如式（6）所示的非線性權：

通過評估每個小模板的光滑程度，調整每個子模板的權重。對于含間斷的子模板，權重自動設置為一個極小的值，從而消除不光滑模板的貢獻。對于5 階WCNS-E 格式，具體計算方法如下所述。

將插值變量到特征空間：

對特征進行高階插值：

式中：f、S分別為一階導數及二階導數的近似。

ω為權重，定義為：

式中：ε是為避免分母為零而加的小量，ε=10－6。ISk為模板的光滑性度量，由式（12）計算：

CLk和CRk為優化權重：

為了匹配高精度的空間離散格式，時間離散通常采用顯式Runge-Kutta 格式。文中采用3 階精度的TVD Runge-Kutta 格式，如式（14）所示：

2 驗證算例

2.1 Blast wave 測試

Blast wave 測試[20]描述的是兩道相向運動的強間斷相互作用，是考核程序的魯棒性和間斷捕捉能力的經典算例。該算例的初始條件設置如式（15）所示：

計算域為x∈[0,1]，在x=0 與x=1 位置施加反射邊界條件。計算域采用均勻網格離散，網格點數為N=400，計算無量綱時間t=0.038?！熬_解”為WENO5-JS 格式，用網格點N=2500 計算獲得。

測試結果（如圖1 所示）表明，采用不同的WCNS格式，可以看出當前程序求解的結果都能穩定且準確地捕捉強間斷，并且與參考解吻合良好。證明了程序采用數值格式的在處理強間斷能保持很好的魯棒性。

2.2 NACA0012 翼型繞流算例

該算例旨在考察程序模擬復雜幾何構型流場的能力。自由來流條件設置如下：來流馬赫數為0.8，基于弦長的雷諾數為500，攻角α=10°。NACA0012翼型弦長為1，周向分布377 個網格點，徑向74 個網格點，近壁處最小法向網格距離為1.7 mm。翼型前緣網格區域擴展到18 倍弦長，后緣到35 倍弦長。出口和入口都采用遠場黎曼邊界條件，翼型表面采用無滑移、無穿透壁面邊界條件，展向依然采用滑移壁面邊界條件，翼型區域總網格大小為377×74×21，后緣到網格最右端網格大小為68×74×21。計算結果如圖2 和圖3 所示。

圖1 Blast wave 測試算例Fig.1 Blast wave case

圖2 NACA0012 翼型馬赫數分布Fig.2 Mach filed around NACA0012

圖3 不同求解器得到的NACA0012 翼型表面壓力曲線分布對比Fig.3 Comparison of the curve of pressure on the surface of NACA0012 by three different solvers

圖2 為HFS 計算得到的機翼流場馬赫數等值線分布，可以看出，流場中的等值線非常光滑平順。圖3 顯示了機翼表面的壓力系數，并將HFS 的計算結果與NASA 的OVERFLOW 程序的計算結果和實驗結果進行了對比。結果顯示，HFS 的結果與OVERFLOW和實驗測量結果均符合良好，驗證了HFS 程序計算可壓縮黏性流場的能力和準確性。

2.3 鈍錐表面脈動壓力分析

該算例旨在考察模擬鈍錐邊界層流動的能力。自由來流條件設置如下：來流馬赫數為6，單位長度雷諾數為107/m，頭部半徑為1 mm，壁面擾動邊界條件見李新亮等[21]，擾動幅值為0.005。沿鈍錐母線方向分布3600 個網格點，周向分布1500 個網格點，沿壁面法向分布360 個網格點，在鈍錐尾部部署200 個拉伸網格用于吸收人工截斷帶來的數值擾動，壁面采用無滑移+吹吸氣擾動邊界條件，遠離壁面采用無反射邊界條件。從圖4 和圖5 中可以看出，鈍錐外部流動分為3 個明顯的區域，初始為層流階段，依次經歷轉捩和完全發展湍流，后續將繼續利用HFS 求解器研究鈍錐擾流的脈動壓力的時空分布規律。

圖4 鈍錐表面脈動壓力云圖Fig.4 Fluctuation pressure on the surface of a blunt-body

圖5 鈍錐外部的脈動壓力云圖Fig.5 Fluctuation pressure near the surface of a blunt-body:a) laminar flow stage;b) transition stage;c) turbilent stage

3 結論

1）開發了一套可壓縮流動脈動壓力求解器，并通過幾個算例驗證了數值格式、邊界條件等的正確性。

2）研究了鈍錐擾流，捕獲了鈍錐邊界層的層流、轉捩和完全發展湍流等流動現象，獲得了壁面脈動壓力的空間分布。